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1.  Introduction 

In  the  last  decade  several  very  fast  methods  have  been 
developed  for  the  direct  solution  of  the  quite  special  systems  of 
linear  algebraic  equations  which  arise  when  Laplace's  or 
Helmholtz's  equation, 

-Au  +  cu  =  f  ,    c  =  constant  , 

is  solved  by  standard  finite  difference  methods  on  certain  simple 
regions  in  the  plane.   The  best  known  of  these  methods  are  due  to 
Hockney  [23,25]  and  Buneman  [5],  see  also  Fischer,  Golub,  Hald, 
Leiva  and  Widlund  [15]  and  Section  6  of  this  paper,  in  which  a 
so-called  Fourier-Toeplitz  method  is  developed.   For  further 
studies  of  such  methods  see  Dorr  [14]  and  Buzbee,  Golub  and 
Nielson  [  8 ] . 

All  these  methods  can  be  regarded  as  efficient  computer 

of  the  separation  of  variables  method.   That  method 
can  only  be  used  for  regions  which,  after  a  possible  change  of 
independent  variables,  are  rectangular  and  for  differential  opera- 
tors of  a  special  form.   Typical  examples  are  Laplace's  and 
Helmholtz's  equations  in  Cartesian  coordinates  on  rectangular 
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2 
of  equations  with  n  unknowns  arising  from  Laplace's  equation  on  a 

square  region  using  n  uniformly  distributed  mesh  points  in  each 

direction,  only  n  (l+(?(l/n)j  storage  locations  are  needed,  while 

o 
the  operation  count  is  Cn  log?n  (l  +  o(l)).   Here  C  is  a  constant 

less  than  ten. 

These  fast  methods  can  also  be  used  for  elliptic  problems  on 
regions,  or  with  boundary  conditions,  which  do  not  allow  for 
separation  of  variables,  provided  that  the  differential  operator 
allows  for  separation  of  variables  on  a  different  region  or  with 
different  boundary  conditions,  see  Buzbee,  Dorr,  George  and  Golub 
[ 7  ],  Buzbee  and  Dorr  [  6  ]  and  Hockney  [26].   This  last  condition 
amounts  to  a  requirement  that  the  differential  operator  can  be 
split  into  the  sum  of  commuting  operators,  each  of  which  contains 
derivatives  with  respect  to  only  one  variable. 

We  will  now  briefly  describe  how  a  problem  can  be  treated  on 
a  bounded  region  O  which  does  not  allow  for  separation  of  vari- 
ables.  We  will  limit  our  discussion  throughout  to  two  independent 
variables.   The  region  O  is  imbedded  in  a  rectangle,  an  infinite 
parallel  strip  or  the  whole  plane.   A  uniform  mesh  is  imposed  on 
the  enlarged  region.   An  expanded  linear  system  of  algebraic  equa- 
tions is  derived  which  has  a  reducible  matrix,  see  Section  J. 
This  new  matrix  contains  the  matrix  of  our  original  problem  as  an 
irreducible  component.   The  resulting  matrix  is  a  rank  p  modifica- 
tion of  a  problem  which  allows  for  the  use  of  a  fast  solver.   Here 
p  denotes  the  number  of  boundary  mesh  points  of  the  original 
finite  difference  equation.   For  two-dimensional  problems  we  thus 
have  a  value  of  p  of  the  order  n.   The  problem  can  be  solved  with 


the  aid  of  the  Woodbury  formula  or  one  of  its  variants,  see  Buzbee, 
Dorr,  George  and  Golub  [7  ].   Similar  ideas  are  also  discussed  in 
George  [18],  Hockney  [26]  and  Martin  [35]  •   It  is  interesting  to 
note  that  Hockney  and  Martin  present  their  algorithms  in  terms  of 
capacitance  matrices,  an  idea  borrowed  from  potential  theory, 
rather  than  in  terms  of  the  Woodbury  formula,  which  can  be 
regarded  merely  as  an  algebraic  identity.   We  note  that  Buckner, 
see  Todd  [^5],  has  pointed  out  a  relation  between  potential  theory 
and  the  method  of  tearing  which  is  closely  related  to  the  Woodbury 
formula.   For  a  further  discussion  of  previous  work  see  Section  8. 

The  work  reported  in  this  paper  on  the  solution  of  the 
interior  Dirichlet  and  Neumann  problems  grew  out  of  an  observation 
of  a  formal  analogy  between  the  Woodbury  formula  and  a  classical 
solution  formula  for  the  Neumann  problem  for  Laplace's  equation  as 
presented  in  Courant-Hilbert  [12],  Garabedian  [17]  and  Petrowsky 
[40] .   In  this  potential  theoretical  approach,  which  goes  back  to 
Neumann  and  Fredholm,  an  Ansatz  is  made  in  terms  of  a  single 
layer  potential.   The  charge  density  is  then  found  by  solving  a 
Fredholm  integral  equation  of  the  second  kind.   This  operator  has 
a  simple  zero  eigenvalue  but  it  is  bounded  and  has  a  bounded 
inverse  on  a  subspace  of  codimension  one.   The  close  analogy  be- 
tween the  integral  operator  and  the  capacitance  matrix  which 
app  i  ^dbury  f       sugge  se 
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practice.   We  have  therefore  chosen  to  solve  the  pxp  linear 
system  of  equations  by  the  conjugate  gradient  method.   The  motiva- 
tion for  the  choice  of  this  method  in  this  context  is  detailed  in 
Section  J.      The  use  of  this  method  results  in  considerable 
savings  compared  to  previous  implementations  of  the  capacitance 
matrix  method  in  cases  where  the  number  of  variables  is  large  and 
only  one  or  a  few  problems  are  solved  for  a  given  region  O.   Our 
interest  in  the  conjugate  gradient  method  began  when  we  searched 
for  an  explanation  for  the  rapid  convergence  reported  by  George 
[18]  with  his  so-called  iterative  imbedding  algorithms,  see  further 
Section  8. 

When  we  turn  to  the  Dirichlet  problem,  we  find  that  the  use 
of  the  Woodbury  formula  corresponds  to  the  attempt  to  solve  the 
original  differential  equation  in  terms  of  a  single  layer  poten- 
tial.  That  approach  is  known  to  give  rise  to  a  Fredholm  integral 
equation  of  the  first  kind  and  thus  an  ill-posed  problem.   While 
it  is  possible  to  prove  that  the  resulting  capacitance  matrix  is 
non-singular  for  any  finite  value  of  p  the  ill-posedness  of  the 
continuous  problem  is  reflected  in  a  growth  of  the  condition 
number  of  the  capacitance  matrices  when  p  increases.   In  our 
experiments,  as  reported  in  Section  9,    we  found  that  the  condition 
numbers  can  become  large  and  also  that  they  vary  rather  erratically 
from  case  to  case.   This  might  lead  to  a  loss  of  accuracy  of  the 
algorithm  as  previously  implemented.   The  conjugate  gradient 
method  also  becomes  quite  uneconomical  in  this  context. 

In  the  continuous  Dirichlet  problem,  the  proper  Ansatz  is  a 
layer  of  dipoles.   Changing  our  Ansatz  to  a  finite  difference 


analog  of  a  dipole  layer,  the  capacitance  matrices,  in  our 
experience,  are  again  quite  well-conditioned  and  the  conjugate 
gradient  method  performs  very  satisfactorily.   We  note  that  this 
successful  Ansatz  falls  outside  the  algebraic  framework  of  Buzbee, 
Dorr,  George  and  Golub  [ 7  ] .   Our  treatment  also  differs  from 
theirs  in  the  Neumann  case  for  Laplace's  equation  in  that  we  allow 
the  capacitance  matrix  to  become  singular. 

A  potential  theory  for  discrete  Poisson  problems  is  developed 
in  Sections  3-5  •   While  our  theory  is  mostly  formal  it  has  served 
very  well  as  a  guide  in  choosing  a  proper  Ansatz  for  the  Dirichlet 
problem  and  certain  scale  factors.   A  rigorous  mathematical 
analysis  of  our  method  will  be  given  in  certain  cases  in  the  thesis 
of  Shieh  [42] . 

Another  problem  in  the  use  of  the  method  as  previously 
developed  is  the  considerable  effort  expended  in  generating  the 
capacitance  matrix.   A  solution  has  been  outlined  briefly  in 
Widlund  [46].   Taken  together  these  improvements  result  in  an 
operation  count  of  Cn  logpn  (l  +  o(l))  compared  to 
Cpn  log?n  (l  +  o(l))  for  the  algorithm  of  Buzbee,  Dorr,  George  and 
Golub  [  7  ] .   The  basic  idea  involved  in  the  fast  generation  of  the 
capacitance  matrix  is  the  use  of  translation  invariance  which  can 
achieved  by  imposing  a  periodicity  cone      as  a  boundary 
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made  independently  by  Peskin  [39]  and  used  in  his  calculations  of 
the  blood  flow  in  the  heart. 

In  our  numerical  experiments,  reported  in  Section  9,    we  have 
used  the  Fourier-Toeplitz  method  for  the  separable  problems.   Other 
fast  Laplace  or  Helmholtz  solvers  could  equally  well  have  been 
chosen.   While  the  experiments  reported  in  this  paper  have  been 
carried  out  using  the  five-point  formula  and  the  well-known 
Shortley-Weller  approximation  to  the  boundary  conditions,  we  be- 
lieve that  our  version  of  the  capacitance  matrix  method  should  work 
equally  well  for  other  difference  schemes.   A  series  of  experiments 
with  a  high  order  method  suggested  by  Kreiss  has  begun.   The  number 
of  conjugate  gradient  iterations  and  the  over-all  performance  of 
the  algorithm  is  virtually  unchanged.   Results  from  these  experi- 
ments and  a  description  of  the  method  will  be  given  in  Pereyra, 
Proskurowski  and  Widlund  [4l] .   It  should  also  be  noted  that  a  fast 
Helmholtz  solver  can  be  used  to  a  great  advantage  in  the  iterative 
solution  of  more  general  elliptic  problems,  see  Concus  and  Golub 
[10]  and  Bartels  and  Daniel  [2].   Similarly  the  use  of  fast 
Helmholtz  solvers  for  transonic  flow  calculations  has  been  advoca- 
ted by  Jameson  [29]  and  Martin  [36,37] .   We  have  also  carried  out 
preliminary  work  to  modify  our  method  to  eigenvalue  problems. 
Dianne  P.  O'Leary  and  the  second  author  are  exploring  the  exten- 
sion of  a  variant  of  our  method  to  three  dimensions. 
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2.  Certain  results  from  classical  potential  theory 

We  give  only  a  brief  review  of  a  few  results  of  potential 
theory.   For  a  detailed  exposition  see  the  references  given  in  the 
previous  section.   We  begin  by  defining  the  potential  ^  resulting 
from  a  charge  distribution  p  on  a  boundary  curve  do, 

Y(x)   =  (l/ir)/  pU)  log  (l/r)dsU)  . 
dO 

We  will  assume  throughout  that  the  boundary  dO  is  sufficiently 

2  2  2 

smooth.   Here  x  =  (x1,x2),  i  =  (i1,i2^    and  r  =(xl~^l^  +(x2~^2^    ' 

We  note  that  (l/2rr)  log  ( 1/r )  is  a  fundamental  solution  of 
Poisson's  equation,  i.e. 

-A(l/2rr)  log  (1/r)  =  6(r)  , 

where  6(r)  is  the  delta  function.   Similarly  the  potential  n    of  a 
dipole  density  m.  on  dO  is  defined  by 

V(x)  =  UA)/  uU)(d/dn^)  log  (l/r)ds(£)  . 
dO 

We  adopt  the  convention  that  the  normal  direction  of  So  is  towards 
the  interior  of  the  region  O  in  which  we  want  to  solve  our  problem. 
We  note  that  the  normal  derivative  in  the  formula  for  °)Y   is  taken 
with  respect  to  the  Greek  variables.   Denote  by  1A  e  limit  of 

?A  when  x  approaches  dO  t     he  interior  and  b\ 
spondinr  identity,  one 

dV7^ 
•'  Vdn  and  .Thus, 


and 


a"2/(+)/an  -   (  +  )    p  +    (1/tr)  J     p(^/Snx)    log    (l/r)ds    , 

ao 

<?|/(  +  )   =   (1)  H  +    (l/ir)J     u(a/an^)    log    (l/r)ds    , 
a?r/an  =   W/bn   . 


ao 


The  Neumann  and  Dirichlet  problems  can  be  reduced  to  Fredholm 
integral  equations.   For  the  interior  Neumann  problem,  we  make  the 
Ansatz 

u(x)  =  (1/2jt)  JJ  f(|)  log  (l/r)d£  +  (1/tt)  J     p(i)    log  (l/r)ds(£) 


O 
us(x)  +  Y(x)    , 


ao 


for  the  solution  of 


-Au  =  f  ,    x  e     O 


du/an  =  g  ,    x  e  ao 


The  first  term,  uc,  is  called  a  space  potential  term.   The  boundary 
condition  is  satisfied  by  choosing  p  such  that 

b1u~/bn   =  -p+  [1/tt)  J     p(a/anxj  log  (l/r)ds 

ao 


(2.1) 


=  g  -  (a/an)uc 


ao 


This  equation  can  be  written  as 
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(I-Kjp  =  g  , 

where  K  is  a  compact  operator  defined  by  the  integral  above.   The 
equation  is  a  Fredholm  integral  equation  of  the  second  kind  and 
thus  a  well-posed  problem.   It  has  a  zero  eigenvalue  and  is 
solvable  if  g  has  a  zero  mean  value.   This  condition  follows  from 
the  standard  Fredholm  theory.   It  can  also  be  understood  as  the 
condition  that  in  steady  state  heat  flow,  the  energy  liberated  by 
heat  sources  in  O  must  be  balanced  by  the  heat  flow  across  the 
boundary. 

If  we  now  attempt  to  use  the  same  single  layer  Ansatz  for 
the  Dirichlet  problem  we  obtain  a  Fredholm  integral  equation  of 
the  first  kind.   It  has  the  form 


(1/tt)  J      p  log  (l/r)ds  =  g  -  u£ 

do 


SO 


and   is   an   ill-posed  problem.      The  kernel   is  now  symmetric,    a  fact 
which  might   initially  appeal   to  a  numerical  analyst. 
The  Ansatz 

u(x)    =    (l/2w)JJf   log    (l/r)d£  +  (1/tt)  J     u(d/dn^)    log    (l/r)ds 
O  do 

=  us(x)   +  W(x)    , 

provides  a  correct  a;  ■  Dirichlet  pre     .   If 

>tes  '  tin 


'  g- 


S 


3. 


This  is  a  well-posed  problem  of  the  form 

(i  +  kMh  =  g 

T 
where  K  is  the  adjoint  of  the  compact  operator  K  introduced  above, 

We  remark  that  the  exterior  Neumann  and  Dirichlet  problems 

similarly  give  rise  to  the  equations 

(I  +  K)p  =  g 
and 

(i-ir)u.  =  g 


respectively.   This  is  important  in  developing  the  existence  and 
uniqueness  theory  for  the  elliptic  problems  by  using  the  Fredholm 
theory  of  compact  operators. 

We  conclude  this  section  by  remarking  that  the  fundamental 
solution  (1/2tt)  log  (1/r)  can  be  replaced  by  any  other  fundamental 
solution  of  the  Laplace  operator.   In  particular,  we  could  have 
used  the  Green's  function  for  a  problem  which  is  periodic  in  one 
or  two  directions.   The  periods  must  then  be  chosen  such  that  the 
region  O  is  a  subset  of  the  interior  of  a  fundamental  region  of 
periodicity,  i.e.  we  must  use  a  wide  enough  infinite  parallel  strip 
or  large  enough  rectangle.   This  observation  is  quite  important 
when  we  implement  our  ideas  computationally. 
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3,  The  imbedding  of  discrete  Poisson  problems 

In  the  next  three  sections  we  will  develop  a  similar,  formal 
potential  theory  for  discrete  Poisson  problems.   For  the  purpose  of 
exposition  only,  we  limit  our  discussion  to  a  case  with  the  same 
uniform  mesh  size  in  the  two  coordinate  directions  and  to  Laplace's 
equation,  i.e.  c  =  0.   The  treatment  of  that  case  is  actually  more 
complicated  than  that  of  c  >  0.   Our  method  can  also  be  extended  to 
negative  values  of  c.   This  requires  a  modification  of  our  current 
fast  Helmholtz  solver  for  the  separable  problem.   We  have  carried 
out  preliminary  numerical  experiments  with  such  problems. 

The  Laplace  operator  is  replaced  by  a  finite  difference 
approximation  such  as  the  five-point  formula.   The  fundamental 
solution  (1/2tt)  log  (1/r),  used  in  Section  2,  will  then  be  replaced 
by  the  discrete  Green's  function  for  the  surface  of  an  infinite 
cylinder,  the  surface  of  a  torus  or  the  whole  plane.  Very  effi- 
cient and  accurate  Poisson  and  Helmholtz  solvers  exist  for  such 

problems,  see  Section  6.   We  will  denote  by  B  the  matrix  repro- 
of 
senting  the  discrete  Laplacian  -h  Ah,  employing  undivided  diff 

ences,  for  the  region  and  boundary  e. 

We  decompose  the  set  of  mesh  points  into  three  dir      sets 

Oh,  bOh   and  (CO)h-   The  set  Oh  is  the  .  r  mesh  points, 

.  each  of  its  men         all  its  *  he 

ope:       .         aininr  >   •  tne 

set      regular  mes:  '.et  (C  the 

r  and 


produces  values  of  a  mesh  function  even  for  the  points  in  (CO).. 
These  values  are  largely  arbitrary,  a  useless  by-product  from  the 
fast  Poisson  solver.   It  is  likewise  necessary  to  provide  some 
largely  arbitrary  extension  of  the  data  to  the  set  (CO),-   We  note 
that  the  formulas  for  the  continuous  problem,  given  in  the  previous 
section,  can  be  interpreted  as  having  resulted  from  the  extension 
by  zero  of  the  right-hand  side  f. 

For  all  regular  mesh  points,  i.e.  those  in  O,  U  (CO),,  we  use 
the  basic  discretization  formula  of  the  Laplace  operator  which  we 
have  already  selected.   That  formula  cannot  be  used  for  the 
irregular  mesh  points,  i.e.  those  in  dOh>  because  we  must  introduce 
an  approximation  of  the  boundary  conditions,  and  we  must  also  assure 
that  our  solution  on  Oh  U  do,  is  independent  of  the  values  of  the 
data  and  the  solution  on  (CO),.   This  will  be  achieved  by  elimina- 
ting from  the  discrete  Laplacian,  centered  at  an  irregular  mesh 
point,  the  values  of  the  solution  at  its  exterior  neighbors.   Thus, 
in  the  Dirichlet  case,  we  will  combine  the  discrete  Laplacian  with 
interpolation  formulas  which  approximate  the  boundary  condition, 
while  in  the  Neumann  case  we  will  use  difference  approximations  to 
the  normal  derivative  for  the  same  purpose.   Examples  with  full 
details  will  be  given  below.   The  resulting  linear  equations 
corresponding  to  the  irregular  mesh  points  should  be  multiplied  by 
appropriate  scale  factors.   Because  the  choice  of  these  factors 
is  quite  important  for  our  theory  and  the  performance  of  the 
conjugate  gradient  method,  it  will  be  discussed  at  some  length 
below. 
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We  will  denote  by  A  the  matrix  corresponding  to  our  differ- 
ence equations  on  the  entire  mesh.   The  equations  and  unknowns  are 
ordered  in  the  same  way  as  the  regularly  structured  problem  given 
by  the  matrix  B.   The  following  representation  of  the  matrix  A  are 
direct  consequences  of  the  fact  that  the  rows  of  A  and  B  are  the 
same  for  all  regular  mesh  points.   For  the  Neumann  problem  we  write 

A  =  B  -  UVT  , 

while  we  use  the  notation 

A  =  B  + UZT 

for  the  Dirichlet  case.   The  matrices  U,  V  and  Z  have  p  columns, 
where  p  is  the  number  of  points  in  the  set  d(")h.   The  matrix  U 
represents  an  extension  operator.   It  maps  any  mesh  function 
defined  only  on  do.  into  a  function  on  all  mesh  points.   It  re- 
tains the  values  of  the  mesh  function  on  dOh  and  makes  the  re- 

T 
maining  values  equal  to  zero.   Its  transpose,  U  ,  is  a  trace  opera- 

T 

tor,  i.e.  U  maps  any  mesh  function  defined  for  all  mesh  points 

T        T 
into  its  restriction  to  dOh.   The  rows  of  V  ,  and  -Z  ,  are  simply 

the  difference  between  the  :  B   and  A.   One  can 


,  „T 


which  I  ws  co' 

have  been  deleted. 

chc 

I  lngs 

• 


PA  = 


(A±1        0 


^  A21   A22y 


The  submatrix  A,-,  is  the  matrix  for  the  linear  system  of  equations 
which  we  set  out  to  solve.   It  is  easy  to  see  from  the  structure 
of  the  matrix  A  that  the  solution  on  Oh U do,  is  independent  of  the 
solution  and  data  on  (CO  J,.   The  second  diagonal  block  App  can  be 
interpreted  as  a  finite  difference  approximation  to  a  Dirichlet 
problem  on  COc   It  is  therefore  easy  to  verify  that  Ap„  is  non- 
singular  if  our  basic  difference  approximation  is  of  positive  type, 
Similarly  any  convergent  difference  approximation  to  the  interior 
Dirichlet  problem  will  give  rise  to  a  nonsingular  A,-,  and  thus  a 

nonsingular  matrix  A.   For  the  Neumann  case  we  assume  that  the  row 

T 
sums  of  A-,-,,  and  V  ,  vanish  and  that  the  matrix  A,,  has  a  simple 

zero  eigenvalue.   It  is  then  easy  to  see  that  the  matrix  A  also 

will  have  a  simple  zero  eigenvalue. 
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4.  A  capacitance  matrix  method  for  the  Dirichlet  case 

We  will  now  describe  our  method  for  solving  the  system  of 
equations 

Au  =  (B  +UZT)u  =  F  . 

Guided  by  the  continuous  analog  we  try  the  Ansatz, 

u  =  GF  +  2GVDu  . 

The  vector  GF  corresponds  to  the  space  potential  part  of  the  solu- 
tion and  satisfies  BGF  =  F.   The  operator  G  thus  plays  the  same 
role  as  the  integral  operator  defined  by  the  fundamental  solution 
of  the  continuous  problem.   Any  constant  belongs  to  the  nullspace 
of  B  and  the  choice  of  G  is  not  unique.   Our  particular 

choice  of  G  corresponds  cloc     3  the  standard  choice  f 
the  .  We  remark  in  particular  tha      an  be 

defined  for  all  .      details  see  Section  6. 

atz  is  a  discrete  analog  of  the 
potential  from  a  dipole  density.  The  vector  u  has  p  components 
and  is  determin  a  system  of  linear  equations  which 

we  will  d      below.  ^.1  and  cc-      scale 

fac*  rs.   The  :lar  rows 

Lrect 

. 

- 


an  approximation  of  an  undivided  normal  difference  of  the  form 
li  (dv/dn)  +o(h).   Here  h  =  h/cos  a  and  a  is  the  angle  between  the 
normal  through  the  irregular  mesh  point  and  the  closest  coordinate 
axis. 

We  now  use  our  Ansatz  and  compute  the  residual  vector, 

(4.1)  Au  -  F  =  (B  +UZT)(GF  +  2GVDp.)  -  F 

=  (2VD +  2UZTGVD)u  +UZTGF  . 

From  the  properties  of  U  and  V,  it  follows  that  the  correct 
difference  equations  are  satisfied  for  all  x  e  Q,.   To  derive  a 

linear  system  of  equations  for  the  vector  u  we  multiply  equation 

T 

(4.1)  by  the  trace  operator  U  .   It  is  easy  to  verify  that 

U  U  =  I  and  also,  for  our  choice  of  V,  that  U  V  =  I  .   Here  I   is 
P  P         P 

the  pxp  identity  matrix.   We  thus  obtain 

(4.2)  (2D  +  2ZTGVD)u  =  -  ZTGF  . 

This  choice  of  u.  will  make  the  residuals  zero  for  x  e   dO,  .   The 
residuals  will  in  general  not  be  equal  to  zero  for  all  x  e    (CO),. 
As  we  saw  in  Section  j5,  this  is  of  no  importance  because  the 
restriction  of  the  solution  to  O  U  dO,  will  still  provide  a  solu- 
tion of  the  original  problem.   The  matrix  on  the  left-hand  side  of 
equation  (4.2)  is  the  capacitance  matrix  C. 

By  analogy  to  the  continuous  case  we  expect  the  capacitance 
matrix  C  to  be  nonsingular  and  in  our  numerical  experiments  we 
have  consistently  found  C  to  be  very  well  conditioned.   When  we 
attempt  to  prove  that  C  is  nonsingular  we  find  that  the  existence 
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of  a  nontrivial  p-vector  <|>  such  that  C<t>  =  0  implies  that  the 
restriction  of  AGVD<J>  to  Oh  U  dCV  must  vanish  identically.   Hence, 
since  A,,  is  nonsingular,  the  mesh  function  GVD<{>  must  be  zero  for 
all  x  e  Oh  U  dfV«   Conversely  if  there  exists  a  nontrivial  <J>  such 
that  GVDtj)  is  identically  zero  on  Oh  U  dOh>  it  follows  from  the  block 
structure  of  PA  that  C<\>  =   U  AGVDsj)  =  0.   We  have  not  been  able  to 
exclude  the  existence  of  such  a  vector  <J>.   In  our  opinion  tools 
from  mathematical  analysis  are  required  to  resolve  this  problem 
completely.   A  proof  that   is  nonsingular  when  B  and  A  are  non- 
singular,  e.g.  when  c  >  0,  is  given  in  Section  8. 

We  would  like  to  be  able  to  interpret  equation  (4.2)  as  a 
discretization  of  the  appropriate  Fredholm  integral  equation  of 
the  second  k'  I.       pie  di         '  n  of  such  an  squati 
would  be  based  on  a  numerical  quadrature  rule  of  the  form 

(4.3)  CK^lDjM^j 

where  r(x,£)  is  ••  solution  of  the 

h.  is  a  local  mesh  1      and  the  values  of  u^  are  approximate 

valuer  .         case  ichlet  data 

-al 

to  choose  a 
sur 

. 
sec 


^, 


length  h. .   The  distance  between  consecutive  irregular  mesh  points 
will  of  course  vary  in  a  highly  irregular  way  and  it  therefore 
seems  appropriate  to  use  a  value  of  h.  which  represents  a  local 
average.   Let  us  consider  a  right  triangle  with  a  slightly  curved 
hypotenuse  made  up  by  a  part  of  the  boundary  of  length  ,/h.   The 
average  distance  between  the  boundary  mesh  points  is  then 
h/cos  a  +  o(h),  where  a  e  [0,tt/4]  is  the  smallest  angle  between  the 
normal  of  do  and  a  coordinate  axis.   We  will  adopt  h/cos  a  as  our 

choice  of  h. . 

1 

T  ** 

Away  from  the  diagonal  the  elements  of  2Z  GV  behave  like 

2h.(d/dru)r  +o(h).   This  follows  from  our  choice  of  h.,  Z  and  V 
and  the  fact  that  divided  differences  of  the  discrete  fundamental 
solution  G  converges  to  the  corresponding  derivatives  of  the 
continuous  fundamental  solution  r(x,  %)    for  x  ^  £,  see  Thomee  [44]. 
Therefore  we  can  simply  adopt  the  choice  D  =  I  .   The  elements  of 
the  capacitance  matrix  C  close  to  the  diagonal  present  particular 
problems.   While  all  off-diagonal  elements  of  the  matrix  resulting 
from  a  numerical  quadrature  are  0(h),  this  is  not  true  in  general 
for  the  elements  of  the  capacitance  matrix  close  to  the  diagonal. 
This  is  a  consequence  of  the  irregular  positions  of  the  points  of 
do,.   Arthur  Shieh  [42]  has  therefore  suggested  that  we  should 
consider  all  elements  of  the  capacitance  matrix  C  corresponding  to 
the  influence  between  pairs  of  points  within  a  distance  of  h  '  ■*    as 
an  entity.   The  resulting  band  matrix  E.  should  then  ideally  be  a 
consistent  approximation  to  the  identity  operator. 

We  will  not  pursue  this  difficult  subject  further  at  this 
time.   We  refer  to  Shieh  [42]  for  a  detailed  theoretical  analysis. 
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However,  we  would  like  to  point  out  that  while  our  description  of 
the  algorithms  so  far  closely  reflects  our  actual  implementation  of 
the  method  we  see  several  alternat :    .   Thus  the  discrete  dipole 
layer  could  be  constructed  differently  by  allocating  a  fraction  of 
the  weights  \i.    to  nearby  discrete  dipoles  located  entirely  out- 
side OhudOh-   This  would  change  the  leading  term  of  the  capaci- 
tance matrix  while  the  elements  far  away  from  the  diagonal  would 
be  virtually  unchanged.   Variations  of  the  elements  of  D  are  also 
possible,  while  retaining  a  local  average  of  one,  without  spoiling 
our  comparison  between  the  equations  (4.2)  and  (2.2). 

In  our  experiments  we  have  always  used  the  five  point  formula 
for  the  regular  mesh  points  and  combined  it  with  the  well-known 
Shortley-Weller  approximation  for  the  irregular  points,  see 
Collatz  [9],  Chapter  5.1  or  Forsythe-Wasow  [ 16] ,  Section  20.9. 
Thus  if  x  £  do.  has  its  eastern  and  southern  neighbors  in  (CO)., 
we  obtain 


(4.4) 


-  (2/(1  +51))uw  +  (2/61+2/62+ch2)uc-  (2/(1 +62))un 


=   h2fc  +  (2/61(l+61))ue  +-  (2/6g(  .       . 


The   sen  Lve   type.  .    =  h,/h  where   h, 

.    along  a  mesh    l  illel  v.  .    of 

1 . 
I 


derive  this  formula  by  using  quadratic  polynomial  extrapolation  to 
find  auxiliary  values  of  u  at  the  exterior  mesh  points  which  are 
next  neighbors  of  x.   These  values  are  then  eliminated  by  using 
the  five  point  formula.   We  note  that  we  can  regard  any  other 
boundary  approximation  from  the  same  point  of  view.   The  sum  of 
the  extrapolation  coefficients  must  always  be  nonzero  and  we  can 
therefore  safely  assume  that  we  will  be  able  to  scale  each 
irregular  row  of  A,  making  each  row  sum  equal  to  one. 

Now  let  c  =  0.   If  we  apply  the  recipe  of  scaling  suggested 
above  we  find  that  equation  (4.4)  should  be  multiplied  by 

(2/51(l+51)  +2/62(l+52))_1  . 

If   only   the    southern  neighbor   is   exterior   the    scale   factor   is 

62(l+62)/2    . 

It  can  easily  be  shown  that  the  diagonal  elements  of  the  resulting 
matrix  A  take  values  in  the  interval  [1,4].   In  many  of  our 
experiments  we  have  used  a  scaling  which  makes  all  the  diagonal 
elements  of  A  equal  to  4.   This  choice  is  close  to  the  recipe 
suggested  by  our  formal  theory  and  it  results  in  a  performance 
which  is  often  slightly  better.   The  use  of  formula  (4.4)  in  its 
unsealed  form  results  in  a  much  poorer  performance,  see  further 
Section  9« 
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5.  A  capacitance  matrix  method  for  ti.  ■  .  -  :jr.ann  case 

In  the  Neumann  case  our  system  of  equations  is  written  as 

Au  =  (B  -  UVT)u  =  F  . 

It  is  solvable  if  and  only  if  the  right-hand  side  F  is  orthogonal  to 
the  left  eigenvector  of  A  which  corresponds  to  the  zero  eigenvalue. 
By  assumption  this  eigenvalue  is  simple.   By  using  the  reducible 
structure  of  A,  described  in  Section  J> ,    it  is  easy  to  see  that  the 
mesh  function  corresponding  to  this  eigenvector  must  vanish  for  all 
x  e  (COL.   The  right-hand  side  F  is  therefore  consistent  regard- 
less of  its  values  on  (CO)h  if  the  data  is  already  consistent  on 

nhuanh. 

Guided  by  the  continuous  analog,  we  try  the  Ansatz 

u  =  GF  +  2GUDp  . 

The  p-vector  p  corresponds  to  the  boundary  charge  distribution  of 
the  continuous  case.   The  extension  operator  U  was  introduced  in 
Section  3  while  the  diagonal      /.  D  contains  certain  scale 
factor  . 

Computing  the  residual  vector  we  obtain 

.1)    Au-F=(B-UVT         Dp)  -F      JD  -2UVTGUD)p  -UV'1  :   . 


which  determines  the  vector  p.   We  note  that  our  formula  is 
virtually  identical  to  the  Woodbury  formula,  see  Householder  [28] - 
The  matrix  of  equation  (5*2)  is  the  capacitance  matrix  C  of 
the  Neumann  problem.   It  is  singular  because  if  it  were  not  our 
procedure  would  provide  a  solution  for  any  data.   The  zero  eigen- 
value is  simple.   Let  us  assume  that  ij),  and  <J>  were  two  linearly 
independent  eigenvectors  of  the  capacitance  matrix  C  corresponding 
to  the  zero  eigenvalue.   Then 


V^UD^  =  DiJk  ,    i  =  1,2 


and  GUD<|>.,  i  =  1,2,  must  be  nonzero  vectors.   These  vectors  are, 
by  a  similar  argument,  linearly  independent.   A  direct  calculation 
shows  that  GUD<j>.,  i  =  1,2,  are  eigenvectors  of  A  with  a  zero 
eigenvalue.   We  have  then  arrived  at  a  contradiction  because,  by 

assumption,  the  zero  eigenvalue  of  A  is  simple. 

T 
The  right-hand  side  V  GF  of  equation  (5«2)  is  consistent  if 

T 
F  is  consistent  for  the  original  problem  Au  =  F.   Let  $     be  the 

left  eigenvector  of  C  corresponding  to  the  zero  eigenvalue.   It 

m        mm 

can  readily  be  shown  that  tJj     =  if)  V  G  is  the  corresponding  eigen- 
vector of  A.   But  ip   F  =  0  implies  <t>T(VTGF)  =  0. 

We  will  show  in  Section  7  that  the  singularity  of  C  will 
cause  little  difficulty  when  we  solve  for  p  by  the  conjugate 
gradient  method.   Any  solution  p  of  equation  (5o2)  will  result  in 
a  zero  residual  for  all  mesh  points  and  thus  a  solution  of  the 
original  problem. 

We  now  attempt  to  choose  the  scaling  of  the  irregular  rows 
of  A  and  the  diagonal  matrix  D  so  that  equation  (5-2)  appears,  as 
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closely  as  possible,  as  a  discretization  of  the  integral  equation 
(5-3)  p  -2  J   (d/dnx)rpds  =  g  . 

Our  arguments  here  are  very  similar  to  those  in  Section  4.   The 
right-hand  side  g  is  by  equation  (2.1)  equal  to  the  normal  deriva- 
tive of  the  space  potential  if  the  Neumann  data  is  equal  to  zero. 

T 
It  follows  from  our  discussion  in  Section  3  that  the  rows  of  V 

are  linear  combinations  of  undivided  difference  approximations  of 

the  normal  derivative  and  the  undivided  discrete  Laplacian.   The 

T 
mesh  function  V  v  is  therefore  essentially  of  the  form 

£h  dv/dn  +  o(h)  where  the  factor  P(x)  is  nonzero.   We  therefore 

scale  the  irregular  rows  of  A  to  make  |3  =  1.   This  choice  makes 

the  right-hand  side  of  (5*2)  a  consistent  approximation  to  that  of 

equation  (5*3)  except  for  an  unimportant  factor  h. 

Adopting  this  choice  of  scaling  of  the  rows  of  A  we  next 

examine  the  elements  of  the  capacitance  matrix  which  are  far  from 

the  diagonal.   They  are  seen  to  behave  like  2h(d/dn  )r  +o(h)  which 

suggests  that  an  appropriate  choice  of  the  elements  of  the  matrix 

D  is  1/cos  a.   Here,  a  e  [0,tt/4]  is  the  angle  between  the  normal 

■ction  and  the  closest  coordinate  axis.   In  our  experiments  we 

have  s      used  as       D  =  . 

As  in  the  Dir       case  we  coul^  r  Ansa-     xn 

to  make  0  ir  discr  is 

■ad  ce  :  ons  of  the  charges  j  mesh 

the  capa 


tried  any  such  procedures  numerically.   An  alternative  Ansatz  of 
this  form  will  lead  to  difficulties  similar  to  those  in  the 
Dirichlet  case.   We  have  thus  been  unable  to  show  that  the  right- 
hand  side  of  the  equation  replacing  (5*2)  is  consistent  except  in 
the  version  of  the  method  which  we  have  used  consistently  in  our 
experiments. 

We  will  now  describe  the  scheme  which  we  used  in  our  computa- 
tions.  Our  point  of  departure  is  the  scheme  (4.4)  for  the 
Dirichlet  problem.   Tne  values  u  and  u  must  be  eliminated  by 
using  some  approximation  to  the  normal  derivative.   Denote  by 
7,  e  [0,7r/2)  the  angle  between  the  normal  direction  at  the  eastern 
point  and  the  x,-axis.   We  have  used  the  formula 

(5-4)   (1-61  tan  7l)uc+61  tan  1±un  -  ug  =  (Bjh/cos  7l)(du/dn)e  . 

Similarly  denoting  by  7?  the  angle  between  the  Xp-axis  and  the 
normal  direction  at  the  southern  point,  we  obtain 

(5-5)       (l-62    tan  72^uc+52    tan  72Uw  "  Us   =    (52h/cos  72)(^u/^n)s    • 
By  combining  formulas    (4.4),    (5.4)    and    (5*5 ),    we   obtain 

-2(1/(1  +61)  +tan  72/(l+62))uw 

+  (2(1  + tan  7l)/(l+61)  +2(1  +  tan  7g  )/( 1  +  6g)  +  ch2  )uc 

-  2(1/(1  +52)  +  tan  7l/(l+61))un 

=  h2fc  -  (2h/((l+61)cos  7l))Ou/Sn)e 

-  (2h/((l  +  52)cos  72))(Su/Sn)s    . 
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When  we  modify  the  formula  in  the  obvious  way  for  the  case  of  only 
one,  a  southern,  exterior  neighbor,  we  obtain 

-(1+2  tan  72/(l  +  62))uw 

+  (2 +2(1  + tan  y2  )/  { 1  +  62  )  +  ch2  )uc  -  (2/(l+62))un  -  ug 
=  h2fc  -2h/(cos  72(l+62))(^u/dn)s  . 

We  note  that  this  scheme  is  of  positive  type. 

According  to  our  recipe  above,  the  scale  factors  for  the 
irregular  rows  of  A  are 

d  =  (2/((l+61)cos  7l)  +2/((l+62)cos  y^)'1 
and 

d  =  cos  72(l+62)/2  , 

respectively.   One  can  easily  show  that  this  results  in  diagonal 
elements  which  vary  be-      >ne  and  three.   Here,  just  as  in  the 
Dirichlet  case  we  ha         ntly  used  a  scaling  which  makes  the 
diagonal  elements  equal  to  four.   This  sea]       thus  quite  close 
to  the  one  suggested  by  our  dis 


6.  The  Fourier-Toeplitz  method  and  the  fast  generation  of  the 
capacitance  matrix 

We  will  now  describe  the  variant  of  the  Fourier-Toeplitz 
method  which  we  have  used  in  our  numerical  experiments  to  solve  the 
Helmholtz  equation  on  the  entire  mesh.   We  will  also  show  how  the 
capacitance  matrix  is  generated  in  our  program  at  a  modest  expense. 
We  have  previously  indicated  that  it  is  important  to  achieve  trans- 
lation invariance.   We  should  therefore  choose  a  region  without  a 
boundary.   In  the  original  paper  on  the  Fourier-Toeplitz  method, 
see  Fischer,  Golub,  Hald,  Leiva  and  Widlund  [15],  the  Dirichlet  and 
the  doubly  periodic  problems  on  rectangles  and  the  Dirichlet 
problem  on  an  infinite  parallel  strip  were  discussed  in  detail. 
Of  these,  the  doubly  periodic  case,  which  amounts  to  solving  the 
problem  on  the  surface  of  a  torus,  would  have  been  an  appropriate 
choice.   However,  we  have  instead  chosen  an  infinite  parallel  strip 
with  a  periodicity  condition,  i.e.  a  problem  on  the  surface  of  an 
infinite  cylinder.   This  variant  of  the  Fourier-Toeplitz  method  is 
somewhat  faster  than  the  doubly  periodic  case  because  an  inter- 
mediary step  employing  the  Woodbury  formula  is  eliminated.   More 
importantly  our  choice  allows  for  a  more  satisfactory  resolution 
of  a  special  problem  associated  with  the  case  c  =  0.   A  third 
alternative  would  be  a  method  of  Hockney's  [25],  p.  178,  which 
allows  for  the  calculation  of  the  solution  of  Helmholtz 's  equation 
on  the  entire  plane.   However,  that  method  requires  four  times  as 
much  storage  and  also  more  arithmetic  operations  and  is  not  there- 
fore appropriate  in  our  present  context.   Hockney's  idea  could 
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however  prove  quite  useful  should  we  attempt  to  compute  the  solu- 
tion of  exterior  problems  without  a  prior  transformation  of  the 
region  into  a  bounded  set. 

The  solution  of  the  five  point  discrete  Helmholtz's  problem 
on  an  infinite  strip  parallel  to  the  x-,-axis  gives  rise  to  a  block 
tridiagonal  linear  system  of  equations  of  infinite  order, 
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represent  the  values  of  the  approximate 
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number  of  indices.   The  matrix  A   is  an  nxn  circulant, 
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For  convenience,  we  assume  that  n  is  even.   The  matrix  A  then  has 

l/2  \T 

the  normalized  eigenvectors  (l/n)    (1,1,  ...,l)   and 

l/2  T 

(l/n)    (1,-1, . . ., -1)   corresponding  to  the  simple  eigenvalues 

p  Q 

2+ch  and  6+ch  respectively,  and  (n-2)/2  double  eigenvalues 
4+ch  -  2  cos  (2rr£/n),  I   =   1,  2,  . . .,  (n-2  )/2,  with  the  eigenvectors 


4" 


,(ij 
H 


<f>l        and   $tt' •      Tnese   are   given  by 


■  U     ; 


^^  =    (2/n)1/2   sin    (k£2rr/n) 


4l,Jk  =    (2/n)1/2    cos    (k^/n)    » 


k  =   0, 1,  .  o . ,  n-1   , 

Z    =   1,2,  ...,  (n-l)/2    . 


A  change  of  basis  corresponding  to  the  diagonalization  of  A 
can  be  carried  out  inexpensively  by  using  the  Fast  Fourier  Trans- 
forms if  n  has  many  prime  factors,  see  Cooley,  Lewis  and  Welch  [11]. 
For  our  program  we  always  assume  that  n  is  a  power  of  2. 

This  change  of  basis,  cf .  Fischer,  Golub,  Hald,  Leiva  and 
Widlund  [ 15] ,  results  in  a  block  tridiagonal  system  with  the  matrix 
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The  diagonal  matrix  D  contains  the  eigenvalues  of  A  .   A  permuta- 
tion of  rows  and  columns  of  this  matrix  which  preserves  symmetry 
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and  groups  the  ,0-th  equation  of  each  block  together  into  one  block 
results  in  a  matrix  which  is  the  direct  sum  of  n  tridiagonal 
Toeplitz  matrices  A  ,  I   =  0, 1, .  ..,n-l,  of  infinite  order.   Here 
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where  \      is  an  eigenvalue  of  A  .   We  will  first  restrict  our 
attention  to  the  case  c  >  0.   All  these  matrices  are  then  positive 
definite.   From  our  assumptions  concerning  the  data  f  we  see  that 
the  corresponding  right-hand  side,  derived  by  using  the  Fast 
Fourier  Transform  on  the  vectors  f^  ,  vanishes  for  j  <  N  and 
j  >  N+  with  N_  and  N+  appropriately  chosen,  we  can  therefore  write 
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tion,  i.e.  requiring 
,  we  find  that 
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By  using  equations  (6.1j  and  (6.2)  for  j  =  N ,  +  1  and  j  =  N_  -  1 
we  obtain 


u  -(i)   -(£)    *(i) 

Vn      UN  +1 
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and 
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We  can  therefore  set  up  a  finite  linear  system  of  equations, 


/• 


^   -1 
-1    A. 


0 


v 


0 


i 

-1 


N 


-l 


N_ 
UN  +1 


U 


V 


N+-l 
(ij 


u 


u 


'ad)  ' 

N 


r 


N  +1 


X-1 
\  n,  y 


This    system   can  be    solved  very   easily   by   an   LU-decomposition 
because    as    is   easily  verified 
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Here  o.  -  \i .  -  \i~.    .   We  note  that  except  for  one  element  in  the 
second  factor  these  bidiagonal  matrices  are  Toeplitz  matrices. 
The  solution  of  the  five  point  discrete  Helmholtz  equation  is  thus 
obtained  by  Fast  Fourier  Transforms  in  the  x2-direction,  by  the 
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while  all  the  others  remain  positive  definite.   Guided  by  the 
continuous  problem,  we  introduce  a  fundamental  solution.   We  choose 
the  function  -|j-k|/2.   We  note  that  we  could  add  any  linear 
function  and  still  obtain  a  fundamental  solution.   With  our  choice 


(6.3) 
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We  note  further  that  no  assumption  of  the  form  £~  "  f, 
needed.   This  is  particularly  important  when  we  generate  a  repre- 
sentation of  the  operator  G  by  using  the  data  f ,  '  =  5..6,  .   The 

k      jg,   km 

formula  (6.3)  can  be  used  to  evaluate  two  consecutive  values  of 
u.  '  and  we  can  then  find  the  remaining  values  by  employing  the 
simple  stable  recursion 

;(0J 


J+2      J+l     J       J+1 


We  now  turn  our  attention  to  the  generation  of  the  capaci- 
tance matrices.   The  matrices  U,  V,  V  and  Z  are  very  sparse  with 
a  fixed  upper  bound  on  the  number  of  nonzero  elements  in  any  row 
and  column.   It  therefore  follows  that  the  number  of  arithmetic 
operations  required  to  generate  a  capacitance  matrix  grows  only  as 
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o 
const,  p  (l+o(l))  if  complete  information  on  the  operator  G  is 

available.   To  obtain  this  information  we  solve  the  Helmholtz 
equation  for  u^  '  '  with  the  data  given  by  f^.   =  5 -o^kO*   In  our 
program  we  take  advantage  of  the  very  special  form  of  the  data 

which  virtually  eliminates  the  need  for  the  first  Fast  Fourier 

( 1  k) 
Transform  step.   If  we  denote  by  uVc"  '  the  solution  resulting 

from  a  single  unit  charge  at  the  point  (j,k)  we  see,  by  the  trans- 
lation invariance,  that 

(j,k)    (0,0) 

u„     =  uv  '  ' 
ly™  i-j,m-k 

We  also  use  the  fact  that  the  solution  u^  '  '    satisfies 

ui°;0)  -  ^'l]    -d    4°;0)  -  A\°l  • 

&,m  -£,m  £,m  £,n-m 


Complete  information  on  G  can  therefore  easily  be  obtained  from 
the  mesh  function  u^  '  '  which  is  generated  by  using  a  simplified 
fast  Helmholtz  solver  only  once. 


7.  The  conjugate  gradient  method 

In  this  section  we  will  briefly  describe  the  conjugate  gra- 
dient method  and  those  of  its  properties  which  make  it  particularly 
attractive  for  solving  the  capacitance  matrix  equations.   We  have 
used  this  method  extensively  as  an  iterative  method  for  the  linear 
systems  of  equations  (4.2)  and  (5.2).   For  a  detailed  theory  of  the 
conjugate  gradient  method  see  Daniel  [13],  Hestenes  [21],  Hestenes 
and  Stiefel  [22],  Kaniel  [31],  Lanczos  [33]  and  Luenberger  [34] . 
In  our  experiments  we  have  used  an  implementation  of  the  method, 
SYMMLQ,  developed  by  Paige  and  Saunders  [38] . 

We  again  denote  the  capacitance  matrix  by  C.   It  is  a  p x p, 
nonsymmetric  and  dense  matrix.   The  conjugate  gradient  method 

requires  a  symmetric  coefficient  matrix  and  we  therefore  work  with 

T 
the  corresponding  normal  equations  Qv  =  b  where  Q  =  C  C.   In  each 

conjugate  gradient  step  we  must  multiply  the  matrix  Q  by  a  vector 

u.   We  never  compute  the  matrix  Q  in  our  program  because  that  would 

o 
require  p  (p+l)/2  multiplications.   Instead  we  first  compute  w  =  Cu 

T 
and  then  Qu  =  C  w.   Each  step  of  the  conjugate  gradient  method 

o 

requires  2p  (l+o(l))  multiplications. 

We  denote  by  v^  an  exact  solution  of  our  system  Qv  =  b  and 
by  v„  the  initial  guess.  In  the  first  step  of  the  algorithm  the 
residual  r„  =  b  -  Qv  is  computed.   This  is  the  steepest  descent 

direction  and  it  is  chosen  as  the  first  search  direction  p-,  .   We 

IT      T 
then  find  the  minimum  v-,  of  the  quadratic  form  ^  v  Qv  -  v  b  along 

the  line  v  =  v  +sp-,.   The  search  direction  p.,  k  _>  1,  is  chosen 

as  that  linear  combination  of  the  residual  r,  -,  =b  -  Qv,  .,  and  the 

k-1     ^  k-1 

previous  search  direction  p,  -, ,  which  makes  the  vectors  p. 

K—  J.  1 
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T 
Q-conjugate,  i.e.   p,Qp.  =  0  for  i  <  k.   The  k-th  approximation  vk 

is  the  minimum  of  the  quadratic  form  along  the  line  v  =  v   ^  +spk« 

No  a  priori  information  about  the  spectrum  of  Q  is  needed.   This 

is  an  advantage  compared  with  the  Chebyshev  and  certain  other 

iterative  methods.   The  iteration  is  terminated  when  the  (>  -norm 

of  the  residual  r,  becomes  less  than  a  certain  tolerance.   For  a 

detailed  description  of  the  termination  criterion  of  the  conjugate 

gradient  method  which  we  have  used  in  our  experiments,  see  Paige 

and  Saunders  [38] . 

It  is  easy  to  verify  that  the  iterates  satisfy 
(7.D  vR  =  v^P^CD^ 

where  P.,  is  a  polynomial  of  degree  k-1.   It  is  known,  see  for 
example  Luenberger  [3^],  that  of  all  iterative  methods  which 
satisfy  the  relation  (7.1)  the  conjugate  gradient  method  is  optimal 
in  the  sense  that 

E(v)  =  \   (v  -v#)TQ(v  -v#) 

is  minimal.   Denote  by  X.,  k^>    i  =  1> !•  the  eigenvalues  and 

eigenvectors  of  Q  and  let 


V0"V» 


nMi 


be  an  eigenvec*  ;»ansion  c  Ltial 

eas 


•*"    w^ko-» 


(7.2)  E(vk)    <  LPk-lf 


for  any  choice  of  a  polynomial  P,  ,  of  degree  k-1.  This  inequality 
remains  valid  if  the  maximum  is  taken  only  over  those  A.  for  which 
a.  ^   0.   The  conjugate  gradient  method  will  therefore,  in  the 
absence  of  round-off,  give  an  exact  solution  in  at  most  q  steps 
where  q  <  p  is  the  number  of  distinct  eigenvalues  of  Q  for  which  a 
coefficient  a.  differs  from  zero.   It  can  be  shown  by  using  the 
inequality  (7.2)  and  a  special  choice  of  the  polynomial  P,  ,  that 

(7.3)   E(vk)  <    (2(1  -  1A)V((1  +  1/A)2K+(1  -  l/y^)2k))2E(v0) 

<  4((1  -  l/vW(l  +  l//^))2kE(vQ)  , 

see  for  example  Daniel  [13].   Here  k   is  the  condition  number  of  Q. 

We  have  seen  in  Section  5  that  equation  (5.2)  has  a  singular 
matrix  when  c  =  0.   It  is  therefore  of  interest  to  consider  the 
performance  of  the  conjugate  gradient  algorithm  in  a  singular  case. 
It  is  easy  to  show  that  the  right-hand  side  b  of  our  normal  equa- 
tions Qv  =  b  is  orthogonal  to  any  eigenvector  <f)  which  corresponds 
to  a  zero  eigenvalue.   Starting  from  a  zero  initial  guess  all  the 
iterates  v,  will  also  be  orthogonal  to  <j).   In  the  absence  of  round- 
off the  entire  computation  will  proceed  in  a  subspace  orthogonal 
to  (j).   The  estimate  (7. 3)  is  then  valid  with  k.   redefined  as  the 
ratio  between  the  largest  and  smallest  strictly  positive  eigen- 
values of  Q.   Our  experience  with  the  conjugate  gradient  method  for 
singular  and  almost  singular  capacitance  matrices  has  been 
excellent,  see  Section  9. 

As  we  pointed  out  in  Section  4,  we  have  not  been  able  to 
prove  that  the  capacitance  matrix  in  equation  (4.2)  is  nonsingular 


36 


when  c  =  0.   As  a  safeguard  we  calculate  the  residual  for  the 
original  nonsymmetric  linear  system  of  equations  (4.2)  and  print  a 
warning  when  its  norm  is  larger  than  expected.   We  have  found  no 
Dirichlet  problems  with  the  constant  c  -  0  for  which  there  are 
indications  that  the  capacitance  matrix  equation  (4.2)  fails  to 
have  a  solution.   This  extra  feature  of  our  program  will  undoubt- 
edly be  useful  when  it  is  extended  to  handle  negative  values  of  c. 
We  note  that  the  solution  of  the  least  squares  problems  with  the 
conjugate  gradient  method  has  been  considered  from  a  theoretical 
point  of  view  by  Kammarer  and  Nashed  [30] . 

In  our  experiments  we  have  found  that  the  estimate  (7.3)  is 
realistic  for  Dirichlet  problems.   For  such  problems  the  capaci- 
tance matrices  are  typically  very  well  conditioned.   We  have  also 
observed  impressive  rates  of  convergence  for  Neumann  problems  with 
very  small  values  of  c.   By  a  continuity  argument  such  problems 
must  have  large  condition  numbers.   To  explain  the  success  of  the 
conjugate  gradient  method  in  these  cases  we  must  consider  the 
entire  spectral  distribution.   The  results  can  be  understood  in 
terms  of  the  estimate  (7.2)  once  we  realize  that  the  small  singular 
values  of  C  are  isolated  and  that  the  singular  values  cluster 
around  a  point  well  removed  from  the  origin. 

We  recall  that  our  choice  of  a  double  layer  Ansatz  for  the 
Dirichlet  problem  was  based  on  the  conje      *hat  1     -.gular 
values  of  the  capacitan  ild  be  d'      ited  similar 

to  those  of  he  classical        al 

>ry.  to  illumina-  a  cas 

which  th 


Kantorovich  and  Krylov  [32],  p.  135,  have  shown  that  for  ellipsoi- 
dal regions  and  Laplace's  equation  the  dipole  Ansatz  for  the 
Dirichlet  problem  leads  to  a  Fredholm  integral  operator  of  the 
second  kind  with  the  singular  values 

1  +7J'   ,     j   =  0,1,  ... 
and 

1  -7J"   ,     J  =  1,2,  ...  . 

Here  -y  =  (a-b)/(a+b)  where  a  and  b  are  the  half  axes  of  the 
ellipse.   The  Neumann  problem,  with  the  single  layer  Ansatz  leads 
to  the  singular  values 


and 


1  -y3      ,  3   =   0,1,  ... 


1  +7J   ,     j  =  1,2,  ...  . 


We  note  that  in  both  cases  there  is  a  very  pronounced  cluster  of 
the  spectrum  at  the  point  one.   Very  rapid  convergence  can  there- 
fore be  anticipated  from  the  estimate  (7»2)  with  an  appropriate 
choice  of  a  polynomial. 

Hayes  [20]  has  shown  that  the  conjugate  gradient  method  con- 
verges superlinearly  for  Fredholm  integral  equations  of  the  second 
kind.   We  have  verified  this  result  by  numerical  experiments  with 
diagonal  coefficient  matrices.   In  our  applications  to  the  capaci- 
tance matrix  method,  we  normally  fail  to  observe  superlinear  con- 
vergence.  The  overall  structure  of  the  spectrum  with  a  pronounced 
cluster  around  the  point  one  is  however  inherited  from  the  continu- 
ous problem,  see  Section  9,    and,  as  we  have  previously  noted,  the 
rate  of  convergence  is  therefore  often  far  better  than  the  estimate 
(7-3). 
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The  single  layer  Ansatz  for  the  continuous  Dirichlet  problem 
on  the  interior  of  an  ellipse  gives  the  singular  values 

7  >  j  =  0, 1,  . . .  . 

A  similar  cluster  of  singular  values  very  close  to  the  origin  can 
be  observed  in  experiments  with  the  analogous  discrete  Ansatz. 
This  explains  the  slow  convergence  and  erratic  behavior  of  the 
conjugate  gradient  method  when  the  Woodbury  formula  is  applied  to 
the  Dirichlet  problem,  see  Section  9« 

In  our  program  we  also  have  an  option  of  using  Gaussian 
elimination  for  equations  (4.2)  and  (5.2)  if  the  capacitance  matrix 
is  nonsingular.   The  factorization  step  of  Gaussian  elimination 
requires  p  /3  (l+o(l))  multiplicative  operations.   The  conjugate 
gradient  method  produces  highly  accurate  solutions  in  a  number  of 
steps  which  remains  virtually  constant  when  p  is  increased, 
conjugate  gradient  method  therefore  essentially  requires  only 
const,  p  operations.   The  constant  depends  on  the  region,  the 
value  of  c  and  the  boundary  condition.   In  our  experience,  its 
variation  is  not  very  large.   Hence  the  conjugate  gradient  method 
should  be  used  for  large  enough  p  if  only  one  or  a  few  sets  of 
data  are  used  for  th  em.   For  exper     al  evidence 

see  Section  9. 


3 .  Previous  work  on  capacitance  matrix  methods 

In  this  section  we  will  examine  some  earlier  work  on  capaci- 
tance matrix  methods  and  also  derive  estimates  of  the  condition 
number  of  certain  capacitance  matrices.   A  method  of  this  type  is 
described  very  briefly  in  a  paper  by  Hockney  [24],  who  credits 
Oscar  Buneman  for  the  idea.   A  detailed  description  of  this  method, 
its  implementation  in  Fortran  and  a  listing  of  a  program  is  given 
in  Hockney  [26] .   The  problems  considered  there  differ  in  certain 
ways  from  those  considered  in  this  paper.   The  five  point  discrete 
Laplacian  is  used  on  a  uniform  mesh  in  a  rectangular  region.   A 
large  variety  of  boundary  conditions,  all  of  separable  type,  are 
allowed  on  the  sides  of  the  rectangle.   A  number  of  electrodes  are 
introduced  in  the  interior  and  on  the  boundary  of  the  rectangle. 
Each  is  represented  by  a  straight  line  segment  to  which  one  or 
several  mesh  points  are  assigned.   Prescribed  average  values  of 
the  potential  on  each  such  set  of  mesh  points  are  obtained  by 
using  a  capacitance  matrix  calculation.   If  the  length  of  a  line 
segment  is  chosen  to  be  zero,  an  electrode  will  consist  of  a  single 
mesh  point.   Although  primarily  designed  for  other  applications, 
Hockney' s  program  can  be  used  for  the  solution  of  an  interior 
Dirichlet  problem  for  Laplace's  equation.   The  capacitance  matrix 
is  then  the  restriction  of  the  solution  operator,  for  the  problem 
of  our  choice  on  the  rectangle,  to  the  set  of  irregular  mesh 
points.   Hockney' s  capacitance  matrices  are  always  symmetric, 
positive  definite,  and  storage  and  time  is  saved  by  using  the 
Cholesky  algorithm.   Hockney 's  method  thus  corresponds  to  a  single- 
layer  Ansatz  for  the  Dirichlet  problem.   In  a  typical  application 
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of  his  program  the  solution  is  required  for  the  entire  rectangle, 
a  fact  which  seems  to  preclude  the  use  of  discrete  dipoles,  cf. 
Section  4.   It  would  not  seem  possible  to  use  his  program  without 
modifications  to  obtain  a  more  than  first  order  accurate  solu- 
tion of  elliptic  finite  difference  equations  on  a  general 
interior  region.   We  doubt  that  the  symmetry  of  the  capacitance 
matrix  can  be  maintained  in  higher  order  cases.   See  discussion 
below. 

The  capacitance  matrices  are  generated  at  the  expense  of  p 
calls  of  the  subroutine  which  gives  the  solution  on  the  rectangle. 
Here  p  is  the  number  of  electrodes.   The  cases  of  Neumann  and 
periodic  boundary  conditions  on  the  rectangle,  which  make  the 
Laplace  operator  singular,  are  discussed  and  the  related  difficulty 
is  well  resolved.   The  capacitance  matrices  are  symmetric  even  in 
these  cases.   By  carefully  examining  his  fast  Poisson  solver, 
Hockney  finds  a  short-cut  for  the  treatment  of  so-called  boundary 
electrodes,  which  are  located  on  two  parallel  sides  of  the 
rectangle.   This  observation  leads  to  the  use  of  two  capacitance 
matrices  for  the  boundary  electrodes  and  the  remaining  electrodes 
respectively.   Once  the  Cholesky  factors  of  the  capacitance  matrix 
for  the  boundary  electrodes  are  computed,  the  correction  due  to 
boundary  electrodes  is  handled  quite  inexpensive      Le  the 

The  report  contains  detailed  inform,-!.        the  accuracy  and 
execution  time  of      rogram  when  IBM  360/91. 

The  reports  b\       .  ;e  an:       [7]  an      ;e 
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results  are  formulated  in  matrix  language.   In  both  papers  the 
imbedding  idea  is  discussed  generally.   Focusing  on  the  first  of 
these  papers,  we  find  that  the  Ansatz 

u  =  B-1F  +B-1UWp 

is  used  when  the  matrix  B  is  nonsingular.   Here  B  is  the  operator 
of  our  choice  for  the  rectangle,  W  a  p  x p  nonsingular  matrix  and 
U  the  extension  operator  introduced  in  Section  3°   A  choice  of 
W  =  I  gives  the  Woodbury  formula.   The  case  of  a  singular  matrix  B 
is  treated  in  a  way  different  from  Hockney's  and  ours.   A  non- 
singular  capacitance  matrix  results  in  that  case  if  and  only  if  A 
is  nonsingular.   An  operation  count  for  this  method  has  already 
been  given  in  Section  1. 

A  symmetric  capacitance  matrix  is  obtained  when  this  single- 
layer  Ansatz  is  used  for  the  Dirichlet  problem  if  no  irregular  mesh 
point  has  a  neighbor  outside  the  closure  of  O  and  W  =  I.   We  see 
no  practical  possibility  of  choosing  a  W  which  will  retain  the 
symmetry  of  the  capacitance  matrix  when  a  more  accurate  boundary 
approximation  is  used.   Symmetry  is  also  lost  for  the  Neumann 
problem  even  for  the  simplest  boundary. 

We  refer  to  Bramble  and  Hubbard  [3]  for  a  proof  that  first 
order  accuracy  is  obtained  if  a  general  boundary  of  a  Dirichlet 

lem  is  approximated  by  a  set  of  mesh  points.   We  also  note 
that,  close  to  the  boundary,  approximations  of  the  gradient  of  the 
solution  obtained  from  such  a  scheme  cannot  be  expected  to  con- 
verge, cf .  Bramble  and  Hubbard  [4] .   This  is  in  contrast  to 
schemes  with  a  more  accurate  treatment  of  the  boundary  condition 
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for  which  asymptotic  error  expansions  exist.   From  these  expansions 
the  second  order  accuracy  of  difference  approximations  to  the 
first  and  at  times  higher  order  derivatives  can  be  obtained,  cf . 
Pereyra,  Proskurowski  and  Widlund  [4l] .   It  should  also  be  noted 
that  the  use  of  isoparametric  finite  element  methods,  cf.  Strang 
and  Fix  [4}],  enjoy  an  increasing  popularity  for  second  order 
elliptic  problems.   In  these  methods  the  boundary  is  approximated 
quite  accurately  by  piecewise  polynomials  of  degree  two  or  higher. 
The  success  of  such  methods  undoubtedly  reflects  the  importance  of 
an  accurate  representation  of  the  boundary  in  many  applications. 

One  can  of  course  also  use  a  nonsymmetric  capacitance  matrix 
in  the  Buzbee,  Dorr,  George  and  Golub  [7]  method.   However,  we  see 
no  reason  to  prefer  their  single  layer  Ansatz  for  the  Dirichlet 
problem  to  our  double  layer  Ansatz.   Only  a  marginally  lower  opera- 
tion count  will  result  from  their  choice  except  in  a  case  where  no 
irregular  mesh  point  has  a  neighbor  in  the  complement  of  the 
closure  of  O  or  an  only  first  order  accurate  approximation  is 
accepted.   In  such  cases  the  savings  in  operations  result  almost 
entirely  from  the  use  of  the  Cholesky  method  instead  of  Gaussian 
elimination  for  the  factorization  of  the  capacitance  mat; 

Buzbee,  Dorr,  George  and  Golub  [7]  also  show  how  capacitance 
matrix  techniques  can  be  used  to  handle  problems  on  unions  of 
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Dorr  [6].   Some  of  these  results  can  be  found  in  Table  8.   The 
main  emphasis  of  their  paper  is  the  application  of  a  capacitance 
matrix  method  to  a  finite  difference  approximation  of  rectangular 
clamped  plate  problems.   Their  method  is  an  improved  version  of  an 
algorithm  derived  by  Golub  [19] •   The  variables  can  be  separated 
for  this  plate  problem  once  the  boundary  conditions  on  two 
opposite  sides  of  the  plate  have  been  suitably  modified.   The 
order  of  the  capacitance  matrix,  p,  is  therefore  twice  the  number 
of  mesh  points  across  the  plate.   To  ensure  flexibility,  the 
matrix  decomposition  method,  cf.  Buzbee,  Golub  and  Nielson  [8],  is 
used  as  a  fast  solver.   This  method  allows  an  arbitrary  number  of 
mesh  points  in  both  directions  but  it  requires  on  the  order  of  n 
operations.   An  efficient  implementation  of  this  method  still 
allows  very  fast  execution.   When  generating  the  capacitance 
matrix,  Buzbee  and  Dorr  take  advantage  of  the  sparsity  of  the  data 
vectors  and  the  fact  that  the  solutions  are  needed  only  at  rela- 
tively few  points.   This  computation  therefore  requires  only  on 

2  "5 

the  order  of  pn  operations  compared  to  pn^  needed  for  a  straight- 
forward implementation  of  the  capacitance  matrix  method  using  the 
matrix  decomposition  method  for  a  general  value  of  n.   We  see  no 
possibility  of  reducing  this  operation  count  substantially  by  a 
particular  choice  of  n  unless  the  problem  is  modified  to  allow  for 
the  use  of  translation  invariance,  cf.  Section  6.   Buzbee  and  Dorr 
also  exploit  obvious  symmetries  of  the  rectangular  plate  problem 
when  computing  the  capacitance  matrix.   In  their  experiments  with 
Laplace's  equation  they  take  advantage  of  a  choice  of  n  =  2  , 
k  integer  by  using  a  faster  Laplace  solver  than  the  matrix  decom- 
position method. 

44 


Once  the  capacitance  matrix  and  its  Cholesky  factors  have 
been  found,  the  solution  of  the  clamped  plate  problem  requires  the 
usual  two  calls  of  the  fast  solver  subroutine  and  the  use  of  the 
backsolving  part  of  the  Cholesky  algorithm  for  a  matrix  of  order  p. 

George's  [18]  paper  contains  a  very  interesting  idea. 
George  notes  that  a  residual  vector  for  the  capacitance  matrix 
equation  Cp  =  b  can  be  calculated  at  the  expense  of  one  call  of 
the  fast  Laplace  solver  subroutine  even  if  the  elements  of  C  have 
not  been  computed.   If  the  capacitance  matrix  is  positive  definite 
symmetric,  the  corresponding  linear  system  can  therefore  be  solved 
by  the  conjugate  gradient  or  the  Davidon-Fletcher-Powell  method. 
George  reports  rapid  convergence  for  these  methods  and  that  the 
number  of  iterations  required  is  proportional  to  /p.   We  note  that, 
in  the  absence  of  rounding  errors,  the  conjugate  gradient  and 
Davidon-Fletcher-Powell  methods  produce  the  same  approximate  solu- 
tions pk  from  the  same  initial  guess  if  the  first  approximation 
of  C   in  the  Davidon-Fletcher-Powell  algorithm  is  the  identity 
matrix.   The  performance  of  the  two  methods  is  also  reported  to  be 
virtually  identical. 

The  Dirichlet  problem  for  Laplace's  equation  is  considered 
in  detail.   It  is  imbedded  in  a  Dirichlet  problem  on  a  rectangular 
region,  a  problem  wl     Ly  str      positive  eigenvalues.   A 
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An  attractive  feature  of  this  method  is  that  the  error  in 
the  boundary  values  is  easily  obtained  in  each  step-   The  error 
elsewhere  is,  by  the  maximum  principle,  bounded  by  the  error  on 
the  boundary.   This  observation  is  used  to  obtain  a  most  conveni- 
ent stopping  criterion. 

George's  numerical  evidence  is  consistent  with  a  conjecture 
that  the  condition  number  of  C  grows  linearly  with  p.   We  will 
prove  a  weaker  result  below.   We  have  computed  the  condition 
numbers  of  certain  families  of  capacitance  matrices  obtained  from 
a  single  layer  Ansatz  for  Dirichlet  problems  and  have  found  a 
linear  growth  with  p,  see  Section  9«   In  that  section  we  also 
report  on  cases  where  the  single  layer  Ansatz  leads  to  an  excessive 

number  of  conjugate  gradient  iterations.   It  should  be  noted  that 

T      T 
we  have  consistently  worked  with  the  normal  equations  C  Cp  =  C  b. 

It  would  be  interesting  to  know  how  the  conjugate  gradient  method 

performs  when  used  directly  for  the  nonsymmetric  linear  system  of 

equations  obtained  by  using  a  single  layer  Ansatz  and  the  Shortley- 

Weller  approximation  for  a  general  region. 

We  will  now  give  an  estimate  of  k(C),  the  condition  number 

of  C  induced  by  the  Euclidean  vector  norm.   We  will  only  treat 

cases  in  which  both  A  and  B  are  nonsingular.   If  we  ignore  the 

T  -1 
factors  2  and  D,  we  can  write  the  capacitance  matrix  C  as  U  AB  U 

for  a  single  layer  Ansatz.   We  need  an  upper  bound  for 
max  |  Cx  |  2  =  max  |  U  AB~  Ux  |  2  ,    |xL  =  1 

and  a  lower  bound  for 

min  |CxL  =  min  |u  AB  Ux  L  ,    |xL  =  1  . 
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Such  bounds  are  provided  by 


and 


max  |UTAB-Iyl2  »    |y l2  =  1 


min  |UTAB_1y|2  ,    |y|2  =  1  , 


respectively  which  in  turn  can  be  estimated  by 


and 


Thus 


IaL-Ib-1!. 


1/(|A^|2|B|2) 


k(C)  <  kU)k(B) 


For  the  method  considered  by  George,  both  A  and  B  have  condition 
numbers  on  the  order  of  h~  .   The  upper  bound  for  |c|2  can  be 
improved  by  the  observation  that  in  this  case,  hC  is  a  quadrature- 
like discretization  of  the  Fredholm  integral  equation  of  the  first 

kind  in  Section  2.   By  using  the  techniques  of  Shieh  [^2],  we  can 

-1 


prove  that  |c|2  is 

«(C)  will  therefor* 

Estimates  of 


and  the  condition  number 
than  const. h~^  or 
;  enable  us  to  distinguish  be- 
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tween  a  Neumann  problem  with  c  =  1,  for  which  the  capacitance 
matrices  in  our  experience  are  uniformly  well  conditi    .  and  a 
single  layer  Ansatz  for  a  Dirichlet  problem.  We  will  now  show,  by 
an  example  provided  by 
imj  ■ 

Consider  two  tridiagonal  mai ■ 


and 
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We  note  that  A  is  a  rank  one  perturbation  of  B  which  results  when 
a  discrete  Neumann  condition  at  one  end  of  an  interval  is  replaced 
by  a  discrete  Dirichlet  condition.   The  order  of  the  capacitance 
matrix  is  therefore  one  and  the  condition  number  equals  one. 
The  inverse  of  B  can  easily  be  computed  because 


where 


L  = 


B 


LL 


N 


-1    1, 


A  straightforward  computation  shows  that 
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/n        (n-1)      ... 
AB"1  =   I  +    (  ^  =   I  +uvT  , 

T  T 

where   u     =    (1,0    ...    0)    and  v     =    (n, n-1,  . . .,  1 ) .      To   find   the    singu- 
lar values,    we   consider 


(I  +uvT)(l  +  vuT)   =  I  +w7j     jW 


T 
where   W  is   a  matrix  with   the    two   columns   v  +(v  v/2)u  and   u.      All 

singular  values  of  AB~   except  two  are  thus  equal  to  one  and  we 

have  only  to  find  the  remaining  two.   These  coincide  with  the 

eigenvalues  of 

1   0 

+  (v  +(vTv/2)u,  u)T(u,v  +(vTv/2)u)  . 
\0 

A  direct  calculation  reveals  that  the  smallest  eigenvalue  is  on  the 
order  of  l/n  while  the  largest  is  on  the  order  of  vr  .   The  condi- 
tion number  of  AB~   therefore  grows  quadratically  with  n.   This 
example  shows  that  the  estimate  of  the  condition  number  of  the 
restriction  of  the  operator  AB~   to  a  subspace  in  terms  of  the 
condition  number  of  AB~"  can  be  v     rude. 

In  order  tc        nat  the  capacitance  mat  ■  iuced 

in  Section  h   are  no:  0,  we  it 

D.     is  ear 
V 
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matrix.   By  using  the  same  arguments  as  in  the  single  layer  case, 
we  find  that 

k(C)  <  k(A)k(B)k(L) 

and  that  C  thus  is  nonsingular.   This  argument  is  again  quite  crude 
and  fails  to  reveal  the  great  power  of  the  double  layer  Ansatz. 

An  interesting  application  of  a  capacitance  matrix  method  to 
the  calculation  of  the  potential  flow  of  an  incompressible  fluid 
past  a  two  dimensional  airfoil  is  given  in  Martin  [35]  .   The 
stream  function  satisfies  a  zero  Dirichlet  condition  on  the  airfoil 
and  is  expressed,  far  from  the  airfoil,  in  terms  of  the  free  stream 
velocity  and  an  unknown  parameter  r,  the  circulation.   An 
additional  condition,  the  Kutta  condition,  allows  the  determina- 
tion of  r-      The  Dirichlet  condition  is  represented  by  an  inter- 
polation formula  which  does  not  seem  to  be  directly  related  to 
classical  finite  difference  techniques.   The  values  of  r  and  the 
appropriate  point  changes  at  certain  mesh  points  inside  and  on  the 
airfoil  are  computed  by  solving  a  capacitance  matrix  equation. 

We  also  note  that  Angel  and  Bellman  [1]  have  explored  an 
imbedding  idea  in  their  work  on  the  numerical  solution  of  elliptic 
problems.   A  matrix  for  a  linear  system  of  equations,  very  similar 
to  our  matrix  A,  is  obtained  by  imbedding  the  original  problem  in 
a  rectangle.   This  system  is  then  solved  by  a  block  Gaussian 
elimination  method.   The  imbedding  simplifies  the  programming 
because  all  blocks  of  the  block  trigiagonal  system  of  equations 
will  be  of  the  same  order. 
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9>    Numerical  experiments 

In  this  section  we  will  report  on  results  from  a  series  of 
numerical  experiments  which  were  carried  out  mainly  on  the  IBM 
360/75  computer  at  the  Royal  Institute  of  Technology  in  Stockholm. 
In  our  experiments  we  have  used  the  program  listed  in  this  report 
as  well  as  some  slightly  different  earlier  versions.   Our  program 
closely  reflects  the  ideas  presented  in  Sections  1-7.   We  report 
on  runs  for  different  values  of  the  mesh  size  and  the  constant  c, 
for  Dirichlet  and  Neumann  boundary  conditions  and  for  three 
different  regions.   Our  main  purpose  throughout  has  been  to  study 
the  performance  of  the  capacitance  matrix  method  as  a  highly 
specialized  linear  equation  solver.   We  have  therefore  worked 
extensively  with  problems  with  no  discretization  error  such  as 
those  whose  solutions  are  polynomials  of  second  degree.   This 
simplifies  the  study  of  the  error  originating  from  the  linear  equa- 
tion solver.   We  are  confident  after  a  series  of  experiments  that 
the  performance  of  the  method  is  virtually  independent  of  the 
character  of  the  data,  the  method  of  extending  the  data  to  (CO)., 
the  width  of  the  strip  and  the  values  of  the  parameters  N  and  N  . 
Some  of  these  experiments  are  further  discussed  bel 

The  first  series  of  experiments  wu  -led  out  to  test  the 
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experiments  for  different  regions  also  showed  similar  pronounced 
clusters  of  the  singular  values  of  the  capacitance  matrices.   We 
note  that  the  capacitance  matrices  in  Table  1  which  resulted  from 
the  single  layer  Ansatz  for  the  Dirichlet  problem  were  fairly  well 
conditioned.   As  expected  they  have  clusters  of  singular  values 
close  to  the  origin.   Other  cases  reported  in  Table  2  are  much  more 
ill  conditioned.   The  erratic  behavior  of  the  spectrum  of  the 
capacitance  matrix  when  a  single  layer  Ansatz  is  used  for  the 
Dirichlet  problem  is  clearly  seen  when  we  compare  the  case  of  c  =0 
and  that  of  c  =10"  .   See  further  the  discussion  in  Section  7.   In 
an  earlier  experiment  the  use  of  the  method  of  Buzbee,  Dorr, 
George  and  Golub  [7]  to  handle  the  case  c  =  0  resulted  in  a  capaci- 
tance matrix  which  had  a  condition  number  l6l  for  p  =  64  and  a 
circular  region.   Thus,  the  use  of  their  method  would  not  appear 
to  resolve  the  difficulties  inherent  in  the  use  of  the  single  layer 
Ansatz  for  the  Dirichlet  problem. 

In  Tables  2  and  3  we  report  on  a  series  of  experiments 
designed  to  assess  the  efficiency  of  the  conjugate  gradient  method. 
The  tolerance  for  the  conjugate  gradient  method  was  set  equal  to 

r 

10"  .   We  note  that  in  most  applications  we  would  be  satisfied  with 
a  cruder  tolerance. 

The  results  reported  in  Table  4  indicate  that  the  norm  of 
the  residual  at  the  termination  of  the  iteration  gives  a  reliable 
estimate  of  the  resulting  error. 

In  Table  5  we  report  on  an  experiment  which  shows  the 
importance  of  a  proper  scaling  of  the  capacitance  matrix.   The 
choice  of  the  scaling  recommended  in  Section  4  resulted  in  19 
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iterations  for  the  case  c  =  0.   We  note  that  the  poorly  scaled 
capacitance  matrices  are  fairly  well  conditioned.   The  large  number 
of  iterations  undoubtedly  reflects  an  unfortunate  distribution  of 
the  singular  values. 

Experiments  were  carried  out  for  the  domain  of  Figure  2  to 
study  the  influence  of  the  data  on  the  performance  of  the  algo- 
rithm.  Dirichlet  problems  with  the  constant  c  =  1  were  used  and 
the  tolerance  for  the  conjugate  gradient  method  was  10~  .   For  a 
problem  with  a  second  degree  polynomial  solution,  12  iterations 
were  required  for  p  =  27  and  15  iterations  were  needed  for  p  =  117. 
The  same  problems  with  random  data  required  13  and  16  iterations 
respectively. 

The  classical  Fredholm  theory  is  directly  applicable  only  to 
regions  with  smooth  enough  boundaries.   Experiments  with  square 
regions,  with  the  sides  of  the  square  at  45  angle  to  the  coordi- 
nate axes  seem  to  indicate  that  our  program  will  not  require  a  more 
than  normal  number  of  iterations  for  domains  with  corners.   We  note 
however  that  such  problems  often  fail  to  have  sufficiently  smooth 
solutions.   This  can  make  the  discretization  error  of  the  finite 
difference  approximation  very  large. 

In  Diagram  1  we  compare  the  rate  of  convergence  predicted  by 
the  first  part  of  inequality  (7.3)  and  that  actually  observed  in 
experiment  . 

In  Table  6         t  the  CPU-time  for  our  program  when  run 
on  an  IBM  360/75  using  a   FORTRAN  H  lev  g  compil   . 

Variations  of  measurements  c       1C"  ;.   We 

.m  recorded 


The  required  storage,  excluding  the  program  itself,  is 
(p+const.)p  +  2n  in  our  implementation.   The  number  of  irregular 
points  p  grows  linearly  with  n.   In  the  case  reported  in  Table  6 
p  is  approximately  equal  to  2n.   We  note  that  the  Helmholtz  solver 

contributes  less  than  20$  to  the  total  execution  time  in  these 

2 
cases.   We  therefore  ignore  the  term  proportional  to  n  log?  n  in 

the  estimate  of  the  execution  time  and  conjecture  that  the  time 

2 
grows  as  const,  p   if  the  conjugate  gradient  variant  is  used. 

When  we  test  this  conjecture  using  the  two  right  most  columns  of 

Table  6  we  find  that  the  ratio  of  the  squares  of  the  values  of  p 

equals  4.36  while  the  execution  time  grows  by  a  factor  of  4.20  for 

the  Neumann  problem  and  by  a  factor  of  4. 36  for  the  Dirichlet 

problem. 

In  order  to  assess  the  efficiency  of  our  program  we  compare 
the  CPU-time  of  our  Helmholtz  solver  with  those  of  other  solvers 
on  a  128x128  rectangular  mesh  as  they  are  given  in  Hockney  [27]. 
The  results  are  given  in  Table  J.      All  times  refer  to  runs  on  the 
IBM  360/75  computer  using  an  H- level  FORTRAN  compiler.   We  note 
that  our  program  is  slower  than  the  others  but  that  this  is  of 
little  importance  in  the  present  context. 

This  and  other  parts  of  our  code  could  undoubtedly  be  made 
faster.   The  program  listed  is  a  modification  of  a  program 
developed  primarily  to  test  the  validity  of  the  discrete  potential 
theory;  relatively  little  attention  was  paid  to  its  speed  of 
execution. 

A  comparison  is  made  in  Table  8  between  the  running  times  of 
two  versions  of  our  program  and  those  reported  in  Buzbee  and  Dorr 
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[6] .   The  runs  with  an  older  version  of  our  program  on  the  CDC  7600 
were  kindly  carried  out  for  us  by  Dr.  Victor  Pereyra  while  he  was 
visiting  the  Lawrence  Berkeley  Laboratory.   When  the  parameters  n 
and  p  were  doubled  the  execution  time  of  our  program  grew  only  by 
a  factor  between  3.52  and  4.10  while  those  of  Buzbee  and  Dorr  grew 
by  a  factor  7-26  or  roughly  as  the  second  and  third  powers  of  p 
respectively.   This  is  as  predicted  by  the  operation  counts.   We 
therefore  expect  our  program  to  be  even  more  competitive  for  larger 
problems.   We  also  note  that  the  ratio  p/n  is  only  one  in  these 
experiments  while  a  more  typical  value  would  be  about  three.   In 
such  cases  with  relatively  few  boundary  points  the  Gaussian 
elimination  option  is  just  about  as  efficient  as  the  conjugate 
gradient  option.   The  experiments  reported  by  Buzbee  and  Dorr  are 
for  quite  simple  geometries  which  would  permit  a  simplification  of 
our  program. 


10.  Outline  of  the  program 

The  program  produces  the  numerical  solution  of  the  interior 
Helmholtz  equation  in  two  variables  with  the  Dirichlet  and 
Neumann  boundary  conditions  by  the  finite  difference  methods 
described  in  Sections  4  and  5.   It  is  written  in  FORTRAN  using 
double  precision  (about  16  decimal  digits  on  the  IBM  36o/75)«   In 
our  experience  it  gives  the  solution  of  the  corresponding  linear 
system  of  algebraic  equations  with  an  accuracy  of  ten  to  fifteen 
decimal  digits  if  the  Gaussian  elimination  option  is  used  or  if 
the  tolerance  in  the  conjugate  gradient  method  is  chosen  small 
enough. 

The  program  consists  of  a  main  program  and  a  number  of  sub- 
routines.  The  subroutines  are: 

INOUT  which  handles  all  input  and  the  output  of  direct  interest  to 
the  user; 

BOUNDR  which  computes  the  required  matrix  elements  etc.  by  using 
the  information  on  the  region  supplied  by  the  user; 
SOLVE  which  is  the  Fourier-Toeplitz  subroutine  for  solving  the 
Helmholtz  equation  on  the  infinite  parallel  strip; 
CAPACI  which  produces  the  capacitance  matrix  C; 
and 

DEKEN  which  is  a  main  subroutine  described  further  below. 
In  addition  the  program  contains  the  subroutines  FORT  and  RFORT 
which  are  Fast  Fourier  Transform  subroutines; 

SYMMQL,  DNORM  and  ATIMES  which  are  double  precision  versions  of 
the  Paige- Saunders  conjugate  gradient  program  and  its  related  sub- 
routines; 
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and 

DECOMP  and  EVLOS  which  are  Gaussian  elimination  subroutines. 

The  user  must  supply  the  necessary  information  about  the 
region  and  the  right-hand  side  as  input.   Details  are  specified  in 
the  comments  in  subroutine  INOUT.   The  boundary  of  the  region  is 
described  in  terms  of  the  coordinates  of  the  irregular  mesh  points, 
their  distances  to  the  boundary  and  the  values  of  the  tangent  of 
certain  angles.   The  program  and  the  data  given  after  the  listing 

of  the  program  are  prepared  for  a  Dirichlet  problem  on  the  region 

2    2 
of  Figure  2.   The  exact  solution  is  (x,+x2-l)/4.   Most  informa- 
tion on  arrays,  etc.  is  passed  through  COMMON.   The  storage 
necessary  for  the  arrays  is  specified  in  the  comments  in  the  main 
program. 

The  user  must  also  supply  values  for  a  number  of  parameters 
in  the  subroutine  INOUT.   They  are: 

CON  which  is  the  real  constant  of  the  Helmholtz  equation; 
L0G2N  which  is  an  integer  such  that  N  =  2**L0G2N  is  the  number  of 
mesh  points  across  the  strip; 

IP  which  is  the  number  of  irregular  mesh  points; 
IH  which  is  an  integer  such  that  IH  <  N/2; 
M  which  is  the  number  of  mesh  points  along  the  strip; 
DELTA  which  is  the  machine  toleran   ; 
and 
AC      :h  is  th      ;sted  toler  :  ..rate  , 

The  c      ution  rail  <        '.he  program  from 

the         Iginating  in  ' 

experience,  within  a  the  paramet 

ACC  . 


Values  must  also  be  given  to  four  logical  variables; 
DIR  with  the  value  .TRUE,  for  a  Dirichlet  problem  and  the  value 
•FALSE,  for  a  Neumann  problem; 

SOL  with  the  value  .TRUE,  if  the  solution  of  the  differential  equa- 
tion is  known  and  the  value  .FALSE*  otherwise; 

ELIM  with  the  value  «TRUE.  for  the  Gaussian  elimination  option  and 
•FALSE,  for  the  conjugate  gradient  option; 
and 

INFUN  with  the  value  .TRUE*  if  for  each  call  of  subroutine  DEKEN 
certain  new  data  should  be  read  in.   Its  default  value  is  .FALSE.. 

The  execution  of  the  main  program  proceeds  as  follows: 

1.  Initialization  of  the  parameters  described  above. 

2.  Boundary  information  is  read  in  and  processed.   Values  of 

IBD  which  is  a  two-dimensional  array  containing  the  coordinates  of 
the  irregular  mesh  points; 

HXX  and  HYY  which  contain  the  distances  from  the  irregular  mesh 
points  to  the  boundary; 
and 

TN  which  contains  the  tangent  of  certain  angles 
are  read  in. 
Values  of 

A  which  is  a  two-dimensional  array  containing  the  matrix  coeffi- 
cients for  the  irregular  mesh  points; 

ID  which  is  a  two-dimensional  array  containing  the  information  on 
the  discrete  dipoles; 
and 
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CFX,  CFY  and  CFC  which  are  one-dimensional  arrays  containing 
weight  factors  for  the  right-hand  side 
are  computed  in  the  subroutine  BOUNDR. 

3.  The  subroutine  SOLVE  is  called  using  the  value  'TRUE'  as  its 
calling  parameter.   This  produces  a  representation  of  the  discrete 
fundamental  solution  which  is  stored  in  the  two-dimensional 
array  F. 

4.  The  capacitance  matrix  is  computed  and  stored  in  the  two- 
dimensional  array  C. 

5.  If  the  logical  variable  ELIM  =  -TRUE,  the  LU  factorization 
of  the  capacitance  matrix  is  computed  and  stored  in  the  two- 
dimensional  array  C. 

6.  The  subroutine  DEKEN  is  called  MAXITR  times.   The  default 
value  of  this  integer  is  one. 

The  execution  of  the  subroutine  DEKEN  completes  the  solution 
of  the  problem.   The  first  five  steps  of  the  main  program  need  not 
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4.  A  collection  of  point  changes  or  discrete  dipoles  is  construc- 
ted. 

5.  The  final  solution  is  computed  by  calling  SOLVE  (•FALSE*)-   It 
is  stored  in  the  two-dimensional  array  F. 

6.  If  the  logical  variable  SOL  =  •TRUE*  the  values  of  the  negative 
logarithm  of  the  absolute  value  of  the  error,  rounded  to  integer 
values,  are  printed  as  in  Figure  2.   If  SOL  =  -FALSE-,  no  informa- 
tion is  printed  about  the  solution  by  the  current  program. 
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Figure  2 


The  last  page  of  output  from  a  run  with  our  program, 
Points  with  a  zero  value  represent  exterior  mesh 
points.   The  strictly  positive  numbers  indicate  the 
number  of  correct  decimal  digits,  rounded  to  the 
nearest  integer,  obtained  at  the  particular  mesh 
points. 
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Table  1 


Neumann 

Problem 

Dirichlet 

Problem, 

Dirichlet  Problem 

Double  Layer  Ansatz 

Single  Layer  Ansatz 

X 

P  =  32 

p  =  64 

P  =  32 

p  =  64 

P  =  32   p  =  64 

0   -0.1 

1 

1 

0 

37 

0.1  -0.2 

0 

0 

17 

16 

0.2  -0.3 

0 

0 

7 

4 

0.3  -0.4 

0 

0 

1 

1 

2 

2 

0.4  -0.5 

1 

1 

0 

0 

1 

1 

0.5  -0.6 

0 

0 

1 

1 

1 

2 

0.6  -0.7 

0 

0 

0 

2 

2 

0 

0.7  -0.8 

0 

0 

10 

24 

0 

0 

0.8  -0.9 

1 

7 

13 

24 

0 

0 

0.9  -1.0 

6 

13 

1 

6 

0 

0 

1.0  -1.1 

11 

22 

4 

4 

0 

0 

1.1  -1.2 

10 

18 

0 

0 

0 

0 

1.2  -1.3 

0 

0 

0 

l 

0 

0 

1.3  -1.4 

0 

1 

1 

0 

0 

1 

1.4  -1.5 

2 

1 

0 

0 

0 

0 

1.5  -  1.6 

0 

0 

1 

1 

1.6  -1.7 

0 

0 

0 

1.7  -1.8 

0 

1 

1 

1.8  -1.9 

0 

1.9  -2.0 

0 

2.0  -2.1 

1 

The  distribution  of  the  singular  values  of  capacitance  matrices  of 
order  p.   The  domain  is  a  circle  and  the  constant  c  in  the 
Helmholtz  equation  is  zero. 
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Table  2 


Neumann 

Problem 

Dirichlet  Problem 
Double  Layer 
Ansatz 

Dirichlet  Problem 
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Ansatz 

P  =  32 

p  =  64 

P  =  32   p  =  64 

p  =  32   p  =  64 
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=  0 

I 
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2.9* 

9      11 
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10       22 

12       20 

c 
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K. 
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44o 

i       15 
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17      35 
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K. 
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67 

12 

6o 

10      14 
4.9     4.4 

13       28 
500      860 
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K 
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13 
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9      11 
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10        21 
10       17 
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Table  3 


Neumann  Problem 

Dirichlet  Problem, 
Double  Layer  Ansatz 

P  =  27 

P  =  56 

P  =117 

P  =  27 

p  =  56 

P  =117 

c  =  0 

I 

K 

13 
3-5* 

15 
3.5* 

15 

16 

3-5 

16 

4.0 

18 

c  =  10"6 

I 
K 

26 
3-5* 

31 
3.5* 

38 
-** 

20 
3.5 

19 
4.0 

22 

c  =  10"3 

I 
K 

20 
3.5* 

23 
3.5* 

26 
** 

16 
3.3 

17 
3.8 

20 
** 

c  =  1 

I 
K 

14 
4.8 

15 
3.3 

16 
** 

12 
2.4 

13 
2.5 

15 

The  number  of  conjugate  gradient  iterations  I  and  the  condition 
number  k:(C)  of  capacitance  matrices  for  the  region  displayed  in 
Figure  2  for  the  different  values  of  c  and  the  number  of  irregular 
mesh  points  p. 


This  condition  number  was  computed  disregarding  a  singular 
value  very  close  to  zero,  see  Section  7« 


** 


Not  available. 


70 


Table  4 


Neumann 

Problem 

Dirichlet  Problem, 
Double  Layer  Ansatz 

P  =  27 

P  =  56 

P  =  27 

P  =  56 

P  =  56 

Norm  of 
residual  on 
exit  from 
SYMMLQ 

2.3 *10~7 

4.o*io" 

6.8-10"7 

6.0-10"1 

4.7*10"' 

Maximum 
norm  of 
the  error 

1.2-10"6 

1.4. lo-6 

-7 
6.0.10 

9.6.10"7 

8.4.10"2 

Normalized 
£2-norm  of 
tne  error 

2.8, 10"7 

4.0. 10"7 

i.9«io-7 

2.2. 10"1 

2.0.10"2 

Number  of 
iterations 

14 

15 

12 

13 

3 

A  comparison   of  the  norm  of  the   residua 
'  thm   and   the    over-all   error   for  pri 
Figure   2  and  with  c  =   1.      The  normalizei 
the   £p-norm  divi':  ire   root   i 

ohu* 


the  conju* 


gradient 


1  £p-norm  of  the  error  is 
)f  the  number  of  points  in 
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Table  5 


Scaling  according 
to  formula  (4,4) 

Scaling  making  diagonal 
element  equal  to  4  +  ch2 

Constant  c 

0 

1 

0 

1 

Number  of  iterations 
with  tolerance 

ACCY  =  1CT6 

>  56 

34 

16 

13 

Condition 
number  of  C 

10 

11.5 

4.0 

2.5 

A  comparison  on  the  effect  of  scaling  of  the  capacitance  matrix 
on  the  performance  of  the  conjugate  gradient  method  for  Dirichlet 
problems  on  the  domain  of  Figure  2.   The  number  of  irregular 
points  was  p  =  56. 
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Table  6 


Number  of  irregular  mesh 
points  p 

27 

56 

117 

Total  number  of  mesh 
points 

16  x  16 

32  x  32 

64  x  64 

Time  for  Helmholtz  solver 
SOLVE 

0.04 

0.15 

0.73 

|Time  to  generate  capacitance 
[matrix  in  the  Neumann  case 

0.06 

0.26 

1.15 

In  the  Dirichlet  case 

0.16 

0.70 

3.06 

Factorization  time  for 
'Gaussian  elimination 

0.08 

0.55 

5.11 

Backsolving  time  for 
Gaussian  elimination 

<  0.01 

0.03 

0.12 

Time  to  reach  10"" 
.tolerance  by  using  the 
'.conjugate  gradient  method 

0.26 

0.96 

4.08 

Time  to  generate  the 
problem  and  other  overhead* 

0.18 

0.53 

2.11 

Total  time  using  conjugate 
gradient  method  in  the 
Neumann  case 

0.50 

1.75 

7.34 

In  the  Dirichlet  case 

0.6 

2.19 

Total  time  with  Gaussian 
elimination,  Dirichlet  case 

•  56 

1.81 

10.62 

Total  time  to  solve  an 
additional  problem  using  t 
Gauss!  . 

0. 

.77 

CPU  time  in  s 
on  an  IBM  360 


the 


Figure  2 


Table  7 


POT  1,  Hockney's  method 

2.55 

XYPOIS,  Buneman's  method 

3.4l 

ODDEVN,  A.  George's  implementation 
of  Buneman's  method 

3-18 

SOLVE,  our  Fourier-Toeplitz 
solver 

3.55 

CPU- time  in  seconds  for  Poisson  solvers  on  IBM  360/75  computers, 
H- level  FORTRAN  compilers  are  used  and  the  mesh  is  128x128. 
The  data  on  Hockney's  and  Buneman's  methods  are  obtained  from 
Hockney  [27]  „ 


74 


Table   8 


Neumann  Problem 

Dirichlet  Problem 

n  =p  =  32 

n  =p  =  64 

n  -p  =32 

n  =p  =64 

Our  current  program 
on  an  IBM  360/75. 
Conjugate  gradient 
option  with 

ACCY  =  10 

0.95 

3.36 

1.14 

4.02 

An  older  version 
of  our  program  on 
an  IBM  360/75. 
Conjugate  gradient 
option  with 

ACCY  =  lO-6 

— 

— 

1.72 

6.19 

The  older  version 
on  a  CDC  7600  with 
the  RUN  76  compiler. 
Conjugate  gradient 
option  with 

ACCY  =  10 

With  Gaussian 
elimination  option 

— 

— 

0.348 
0.338 

1.43 
1.45 

Buzbee  and  I 
Fortran  program,  as 
reported  in  [6],  on 
a  CDC  7600 

— 

— 

O.365 

2.65 

CPU 


L*********  ********* 

C*********  THE  PROGRAM  ********* 

C  *********  ********* 

L*«********S****«************************************************************** 

IMPLICIT  KEAL*8(A-H,0-Z ) 

LOGICAL  ClR,SOL,ELIM, INFON 

CGMMON/iJLKl/     C( 128,128) 

CGMMGN/dLK.2/     F(6<»,64),F0N(6'tf6<»),REl64),Sl(64),LG(64> 

COMM0N/4L*  3/     A(l2d,J),X(12a),M12d),CPX<128), 

*  UY{12d),CFC<12d),HXX(12d),HYYll28),P<123),VU128),V2(128), 

*  *K12d),Si.ALES<Ud),TN(l28,2),IdD(12  8,2),IPS(128),ID(12  8,<t) 
CUMHON/6LKW  ACGY  ,LGN ,CH2, DELTA, h, HSU t IP .  I C, J  0  ,LOG2N,N , M , 

*  max  I TR,  I T ,D1k, SOL.EL IM, INFON 
COHMON/AT/PK  126) 

C 

c***  rhjb  ^KuuKAf  PkODOCES  THE  NOMEkIlAL  SGLGTIGN  GF  THE  HELMHGLTZ 

C*«*  EuOATION  (-LAPLACIAN  «-  CONIU  =  Ft  wlTH  THE  OIKICHLET  GR  NEOMANN 

C***  oUUNUARY  CGNDITIGN,  IN  A  oENERAL  BOONDED  TuG  DIMENSIGNAL  REGION. 

C***  HfcRE    CGiM     lb    A    NON-NEGATIVE    CONSTANT.    A    VARIANT    OF     THE    CAPACITANCE 

L***  MATRIX  METhGU,  DESCRIBED  IN  A  1975  ERDA  REPORT  BY  M. PRGSKOROWSK I 

C***  ANO  G.«lOLUhD,  IS  GSEO  TO  SOLVE  THE  LINEAR  SYSTtMS  OF  EUOATIONS 

C***  rfHICH  ARISE  PRCM  FINITE  UlFPERtNCE  APPROXIMATION  OF  A  CLASSICAL 

C***  TYPE,     wh    ^hhl-K     THE    USER     TO    SECTION    10    OF     THAT     PA^ER    FOR    FORThER 

L*«*  ADVISE  ON  THE  USE  CF  THE  PROGRAM.  THE  PROGRAM  *A S  kRITTEN  IN  FORTRAN 

C***  H  Y  h.PKGSKURCwSKl  AND   aTENSIVELY  TESTED  ON  ThE  IBM  360/75  AT  THE 

C***  ROYAL     INSTITUTE    OP     TECHNOLOGY     IN    STOCKHOLM,     SrfEDEN. 

C**<-  THIS  IS  Tot  NOVEMBER  IS  75  VERSION. 

C 

C«<=*  THt  UitK  MGSI  SUPPLY  HiS  C*\  OAIA  CAkOS  IN  GKUER  TO  CHANGE  ThE 

C»*«  REGION  AND  Tht  DAT;.  ,PACE  MOST  dE  RLUOESTEO  IN  COMMON  FOR  THE 

C***  AkkAYS  LlSTtU  BELO«. 

C 

L«**  iHt     Ph  LL-NSIbTS    GF     A    MAIN    PROGRAM    AM    A    NOMBtR    GF     SObROOTINES. 

C*«*  r  til         JdROUTlNES    AR  E  ;     iNOUTtDEKENi     SOLVE,     BGONDR,     AT1M£S    ANO    CAHACI 

C***  WRITTEN    BY     »i.PK(JSKUROriSKI;     FURT     ANO    RFORT    BY    DR.J.COOLEY; 

L*+*  iYMMLtJ    AND    DNORM    dY    OR. L.PAIGE    ANO    DR. M. S AGNDERS    AND    OECOMP    ANO 

C«*»  EVLOS    fhij'-1    FOKSYTHE    ANO    HOLER'S    BOOK. 

C 

(;*♦♦  \     RECTANGULAR     SObStT     Of      A     UNIFORM     MESH    ON     A    VERTICAL     INFINITE     STRIP 

L***  IS    ObtD    IN     I  nr     LL.MPOTATION.     THE     1-DlKtLTlGN     IS    ACROSS,     THE     J-DIRECTION 

C**«  ALl     ■  ■     1HE     ilR  1  P. 

L 

C*»*  DIMENSIONS    Of     ARRAYS    Rfc*.GESTEO    IN    COMMON. 

C*»*  MULTIDIMENSIONAL    ARRAYS:     F    ANO    FON    GF    OROER    N*M;    C    OF    ORDtR     IP*IP; 

C**»  iBOtTN   01        KOI  R     IR*2;     in    CF  t    IP*4;    A    Of    ORDER  1*»*3 

L«*«  .     ■   •  ,     FY f HXX tHYYt X » NV P ( VI tV^t Wit  SCALES    AND    IPS    GF    ORDt  R     IP ; 

C*»*  REtSI    ANO    Lo    Ul     i  ROEK    Ni 

C 

(.*********<><  **t*«**«***«<>****«*«**«t««c*«««*«**«******«*«*********************** 

C 

C*«*  MIS     RAh!         i       IM      R-UN     RKuGP/tM     1  j     tXECUTEO     ONLY     GNU      f     IK     A     wlvEN 

L***         CHI   H  t    Ul     A    REGlUNt    A     IOUNOARY    CUNDITlUNi     A    VALUE    of     CON    AND    A    MESH. 

C**»         I  .  i  .    PARI    01     IHI     PROGRAM    ! 

C 

C       ll)  'i  ETERS   ARI      ilVEN    THEIR    INITIAI     VALUI     , 

i 


CALL  INCUT  I  I ) 
C 
C   (2)     BOUNDARY  INFORMATION  IS  PROCESSED. 

C 

ULL  1\0UT(2» 

LALL     3GUN0R 
C 
C        (3)  THE     UISLKtlt     FUNDAMENTAL     SOLUTION     IS    COMPUTED. 

C 

l  >\  l  L  SULVEI.TRUE.J 

C 

L        K)     IHE  CAPACITANCE  MATRIX  IS  CUMPOTED  AND  STORED  IN  THE  ARRAY  C. 

C 

I...LL  CAP  AC  I 

c 

C*#*    [(-     „fc  HAVE  A  NtUMANN  PkUBlEM  AND  COIN  IS  ZERO  OK  VERY  SMALL 
C***    the  CONJUGATE  GRADIENT  OPT  lU.N  SHOULD  BE  USED. 

IF(  .NOT.l)IR.ANU.DABS(CH2).LT.1D-6)  ELI M=. FALSE. 

L   151  ThE  CAPACITANCE  MATRIX  IS  FACTORED  IF  THE  GAUSSIAN  ELIMINATION 

C***  uPTION  IS  IN  EFFECT. 

C 

IF  (ELIM)  CALL  LfcCGMPISCALES,  IPS,  IP) 

C 

£$$**************************************************************************♦ 
C 

L*s*    COMPUTE  THE  SOLUTIONS  FUK  MAXITK  D1FFERENI  DATA. 
I. 

DO  200  I  T=l,MAXlTf 

CALL  OEKEN 
200    CONTINUE 

STOP 

END 
C 
C 

£***££*$***$**********************************<--******************** 

c 

c 

SuBRuuTlNE  UEK.EN 

IMPLICIT  kEAL*8(A-H,0-Z ) 

LOGICAL  DIR  , SOL, ELI M, INFUN 

CUMMUN/BLK.I/  Cll28,12d) 

CUMM0N/BLK2/  F(64,64),FUN(6<r,G*t),i<E(b<»),SItt>«t),LuI.6<,) 

CurMON/riLKi/  A<128,3),X(128),*(l2o),CFX(128), 

*  CFY(126),CFC(128),HXX(128)  ,HYY(128),P(128),V1(128),V2(128), 

*  Wl(  128) ,  SCALES! 128) , IN (I  28, 2), IBU( 12 8, 2), IPS(128),ID(128, 4) 
COMMON/ BLK4/     ACCY,CUN,Cn2,DELTA,tl,HSv«,IP,  I0,J0,L0G2N,N,M, 

*  "AX ITK, I T.UIR , SOL, EL IM, INFUN 
C 

L***  THIS    MAIN     SUtJKLUIINE    PRODUCES    THt     SOLUTION    OF     THE    HELMHOLTZ    EUOATICN 

C***  rtlTh    A    SPECIFIC    SET    OF    DATA. 

C 

CALL  IN0UTI31 
C 

C***    COMPOTE  THE  LOAD  FUNCTION. 
C 

DO   30  J=1,M 

UU   30  1  =  1,  N 
30     F (  I  , J) =FUN(  I,J| 
C 
L***    CUMPUTt  THE  SPACt  POTENTIAL. 


c 
c 


35 


40 


50 

c 
c 

c 

c** « 
c 


t>2 

«4 


BS 

c 

(_♦♦♦ 

c 

L 


CALL  SGLVbl  .FALSE.) 

COMPUTE  THE  VECTOR  *  rtHlCH  IS  THE  RIGHT  HANO  SIRE  UF  THE  CAPACITANC 
MATRIX  FUUATIUN. 

CALL  INGUTI4) 

UU  3  5  I*lfIP 

K>IBD< It  II 

L-IHDI  1.2) 

X(l)=CFX(I)*Vlll)*CFY(I)*V2(l)«-CFC(I)*FUN<K,L) 

CONTINUE 

SO4.0D0+CH2 

DU  5  0  1  =  1,  IP 

K=IBD( 1,1) 

L=lbO<  I. 21 

HX=HXX ( 1 ) 

HY=hYY ( 1 ) 

KUR=-DSIGN(  l.OOCFX) 

LUR=-OSIGN(  l.OOO.HY ) 

A  1  =  A  (  1  ,  1  ) 

A2  =  A(  I  ,?) 

IF( .NUT.UlK )  A  3  =  A (I, 3) 

T  =  F(Kf  L)*SC*Al*KK-KORf  L)+A2*F<  K.L-LOR) 

It-  (  DAhSlHX  )  .LE.  1  .000)  GUTO  40 

1P(GI« )  A3=A1 

r=T*A3*F(K  +  K()R,L) 

CUNTINOE 

IP ( UABS(HY)  .Lb.  l.ODO)  bUTO  50 

IF(OIR)  A3=A2 

r*T+A3*F(K.L+L0R) 

*(  I  )  =  x(  1  )-T*HSC 

COMPUTE  THE  SCIUTILN  UF  THE  CAPACITANCE  MATRIX  EJUATIUN. 

IF  ( FL  If  )  CL1C  10b 

CUNJUGATE  CKAU1FNT  CPTIUN. 

DO     6  4     l*l,Ir 

S=0.000 

UU     82     J=  I,  I H 

S*S*C( J»I)*H( J) 

.  I  INUb 
I  (  I  l=S 

.r inuL 
•     i    Lt  II 

■  (  I  )  -     .  jUO 

t  '■  (  I  )  «M  (  I  ) 

nil    MPdl 

Imi     UEFAUll     VALUI     01      I  ml     MAX  Imlh  \      II     -  ATI CN  .       I         | NMl  - 

IS    IP    .  !•*•    THl  F    THE    I   &PA<  I  TANG     MATRi     - 

(.ALL     S>YMMLC(IP»X»»»,P.V1 ;,     ELTAtACCVt  IPtlSTOPI 

■     I F    X    I       ••         i     i  : 

.     ••     l     Lt  IP 
DO 


103 


Uo 


107 

c 
c 

108 

c 

c 
iosi 

c 


110 

c 

c 
c 


115 

c 

118 

C 

C*** 

c 


120 
125 
C 

C 

C 

c 


00  103  J=li IP 
b=S+CI  I  , J)*X( J) 
P (  I  )=S-SLAL  ES(  I  ) 
S  =  000 

00  lOo  1  =  1  f  IP 
S  =  S  +  P(  I  )**2 
S=DSQRT( S/OFLUATl  IP  )  ) 
wklTE(6,107)  S 

FORMAT!'  L2-NCKN  LF  RESIUOAL  CF  CAPACITANCE  MATRIX  EQ  .  =  •  ,  1PD  1  5.  5  ) 
LNO  TtSl 

GUTU  109 

CJNTlNOE 
GAUSSIAN  ELIMNATICN  OPTION. 

CALL  EVLOSI X,H, SCALES, IPS, IP) 

CON T INUE 

COMPUTE  A  FINAL  LOAD  FONCTIOK. 

00  110  J=lfM 
00  110  1=1, N 
F( I iJ)=FUN{ I , J)*HSJ 

IF(DIR)  GOTO  118 

SINGLE  LAYEk  PCTENTIAL 

00  115  1  =  1,  IP 
K=iBD(  I,  1) 
L=IBO{  It  21 
FIK,L)  =  F{K,U+X(I) 

CONT  INOE 
GOTO  125 

CONT  INOE 

DOUBLE  LAYER  PCTENTIAL 

DO  120  1  =  1,  IP 

T  =  TN( It  1) 

K-IB0( It L) 

L=IBD{  1,21 

F(Kf L)-F(KfL)+XtI) 

K=IO(  1,  1) 

L=IO(  I  ,2) 

F(K,L  )=F(K,L)-T*X< i  ) 

K=IDI I ,3) 

L=ID(  I  ,4) 

F(Kf L)=F(K,L)+(T-100)*X(  I) 

CJNI INUE 

CONT  INUE 

COMPUTE  THE  FINAL  SOLOTION. 
CALL  SLLVE [ .FALSE.) 

F  IS  NOW  THE  SCLUTICN  TO  THE  HELMhULlZ  EQUATION. 
CALL  INCUTI5) 

80 


c 

c 

c**** 

c 

c 


RETURN 
END 


fc$#4:$*:«:**  <:**$.«.*<.<..*:*<:*******.***********  ************************ 


C 
(.*** 

c*** 
c*** 
c*** 
c*** 

c*** 

C*** 

c 


SJbROUTlNE  bGONCR 

IMPLICIT  KE«L*bl A-h.U-Z) 

LOGICAL  UlP  ,SGl  ,ELi^,  INFUN 

COMMGN/oLiO/  A(12d,3),X1128),M12d),CFX(12d), 

*  CFY(l2b),CFC(128),HXX(12d)  ,HYY(12b),P<l2b),Vl(128),V2<128>, 

*  *l(12b),SCALES<l^H),TN<L28,2),lBD(12  8,2),IPS(12b),l0(12  8,4) 
COMMON/bLK4/  ACCY,CGNftH2,DELTA,h,HSg, IP,  1 0 , JO , LOG2N, N, M , 

*  MAX  ITk,  1  T  ,L>lfc,  SUL.EL  IM,  INFCN 

THIS  SUbRUUriM-  CUMPUTES  ARRAYS  A  ,  I  0  ,CFX  ,  CF  Y  ,  C  FC  . 

A  CdNTAlNS  THE  CLE FF I C I ENT S  OF  THE  DISCPETE  HELMHOLTZ  OPERATOR 

t-Lk     ImE     l^KtC-ULAk     MESH    POINTS. 

ID  CONTAINS  THE  CULmDIMATES  OF  THE  UlSCkETfc  DIPOLES. 

CFXfCFYtCFC  CONTAIN  THE  COEFFICIENTS  FUR  THE  RHS  OF  THE 

APPROXIMATION  AT  THE  bOUNDARY. 

CO^MuN    UStD     (APART     FRO^    THE     AbOVE    ARR AY S  )  :  I P , DEL T A ,CH2 , H, HSU , D Ik  , 


2b 

30 

C 

C  *** 

c 


c 
c**«- 


sc= 

DO 
I-  I 
J=I 
Hi  = 
hj= 
HX  = 
hY  = 
K=- 
L  =  - 

■ 
DY  = 
CX  = 
CY  = 
It  ( 
Mf  I 
P.JK 
CON 


<»DO+C 
200  I 
bD(  IS 
hD(  IS 
HXX(  I 
riY  Y  (  I 
DAbS( 
UAtiS( 
,  I  ,'. 
USIGN 

U»l  M 
O.uUO 
U.oCC 
HX.Ot 
Tt  U, 
M  All1 
I  INUI 


i-1  ,  IP 

,11 

,2) 

b) 

S) 

HI  ) 

MJ  ) 

(  IDO.hI  ) 

(1D0.HJ  ) 

(  l.ODO.hX  ) 

I  l.ODO.f  V  ) 


.lJlLTA..    .    .HY.OE.DELTAI    c-lUj     io 

2b)     I    . 

Lit  I  \    L.h     h\     Luj     l  HA,-,     JtLTA     IN    POINT     NJ.'.Ii) 


if,  U  I-  i     ■  •  i   i      ■  •   I  MAT  ION 

CC»2G0*SC/1200/I   ■  ♦,,!.u/dy«-Ch2) 

A(  1 S,  I  )  =-CC/(  1  .UDU  +  UX  ) 
..  (  IS, 2  )=-CC/(  1  .OUOH  -  > 

.1  .01)0)       '      i  ■  • .  •        D/OX 

It  Ihy.ii-.I  ..>uu)    CY*CC/(DY*1.0C0)/UV 
Ifl      1.    )      .UTO     IOC 

THE    Nl  PI        •  I    '. 

■  ;  .     .  I) 

All       ,  2  ) 
Ai=U.UbU 

1 1  i  m  ■  .  .  i  . :  . 

1      I  M   I    .,  1  ) 

S*S*CA*<HX*T-1*0I      ) 
,21 


. 


I 


/  I   is ,  1  ) 

■  ■ .  i  >l2J 

I  HI       1 1  AGONAL 


CA=CX*hX*H*OSUK  H  l.0D0+T**2) 
UJ  CONTlNOc 

IFIHY.GT.1. 

1  =  TM  IS,,>) 

S=b+C>  *4hY*  r-1  .CCG ) 

Al=A(ib,l)-CY*HY*T 

CY  =  (  Y*HY*H*CSCK1  (  i.  >  ; 

5  0  LUNT  I 

i  ;   IhX.GT.l  .ODv-  ) 

IHhY.GT.L.  ) 

L  *  *  *  i  L  A  L  L      (L    4< 

S=SC/S 

A(  IS.  1  )  = 

Al IS .2  I  =  A2*S 

A  (  I  S  ,  i  >  =  A  3  *  S 

CX=LX*S 

LY=CY*b 

CC=CC*S 

OuTC     1 S  J 
C***  OUObLL     LAYL* 

lO'j  CuMlNUE 

IO(ISfl)=l+K 

.     (  IS,2)=J+L 

If   [HX.GT.HY ) 

1D(  IS,  3)=I+K 

IDl  IS,4)=J 

GGTC    ISC 

CbMINUF 

iJ{   1  i>  ,  3  )  =  1 

IDl  1  S  ,  4  )  -  J  *  l 

CGuT  INUl 

C  Fa  (  I  S  )  =  CX 

CFY(  1S)=CY 

LI-l  lib.)  =.50C*I-S 

CL-NT  KUt 

KE ! Utf  \ 

END 


R  tP^L  SC.NTATION. 


I       I        1  <.  0 


12>) 


190 


2JJ 


L 
C 
C**: 

c 

c 


*********  *********; 


:*********  *  *********  *  * ********* ; 


******** 


c 

C*** 


SUdROUT  IM      SOI  V[  IRGIN1  ) 
I  m  V  L  I  C  1  T     P.LAl<-8(A-H,()-Z) 
LUG  I  GAL    DIk  ,  SlL  ,  lL  i  N,  1  \rUi\i ,  PL  Ii\  T 

COMMON /BLK 2/    F(  64  .64  >  . FuN«  64  ,  64  )  .Pi  I  64)  .  S  I  (  64)  ,  LG  (64  J 
CJMMG'm/uLKW    ACCY,CCN,CH2.  DELTA, I  ,  \  S.. ■  ,  IP,  I  0 ,  JO  .L0G2N,  N.  M, 
VAXITk, I TfOlP,SGL,QL IV, 1M-UN 

I  Ail     bllVEK    Fl;h     In:      TwU    D 1  ML NS 1UNAI     HLLMHCLT2     EUDAHON    FOR     ALL    CLN>=0 
IN     Mm     INF  INI  I l    PARALLEL     STRIP     hITH    PERIODIC    BOUNDARY    CONDITIONS. 
ThL     ARRAY     F,     Ul      CKCEK     \*M,     *HLkL     N  =  2**N.     k  =  3  .  4  .  •  .  C  L  NT  A  I  KS     RIGHT     HAND 
SIDE    r;     IT     IS    LVCKrtKlTTEN    UURiNG    ThL    tXtCUTION    bY     THE     SULOTICN    U. 

'.      ..     ,    .  'CINT     IS     .iPuC.     IF    I     CONSISTS    ONLY    OF     A    POINT    CHARGE 

AT     (  I  C ,  J  0  )  ;     C  T  F  E  R  w  1  S I      IT     IS     .FALSE.. 

CUVV.uN    USED     :     ARRAYS    FiRE.SI;     VARIAbLFS    LCG2N,CH2 .  DELTA  ,N,  MfHSQt 
10.  JO;     iJi^K  IL  r  [iNE    RFOF  I  . 

PkK      *PgTE     s  INE 


82 


c 
c 
c 

c 


35 
40 

C 
50 

C 

c*** 


70 


00 
90 

c 

c*«=* 
c 
c 

c 

L  *<■  « 

c 


c 

c 
c 

c 


I  10 

I 

c 

1  12 


CALL  PFCPT(RE,LIG2N, SI, J, IFERR) 

IKPGIMIGUTO    50 

FUOklER  TkANSFLR^ATION  OF  F 

00  40  1=1, M 

DO  JO  J=1,N 

KC< J)=F< J, I ) 

CALL  RFCRT(RC,LCC2N,Sl,-2, IFERR) 

DO  35  J=1,N 

FIJ,  I  )  =  KE< J ) 

COM  I  NUE 

GOTO  90 

COM  1NUE 
PUIM  CHARGE  GPTICK. 

IMTIALlY  F  AND  RE  ARE  SET  TC  ZERO 

Ou     7C     [  =  l,N 

kL ( 1 )=0.000 
DO    70    J=1,H 

F(  I  ,  J)=OL>0 
RE (  101  =  100 
CALl     kFLkT(RE,LCG2N,S1 ,-2, IFERK) 

60    J  =  1 »  N 
F(  J, JO)  =  RE( J) 
COM  INOE 
PIN  =  U.  JCO*DATAM  1  .ODOJ/DFLGAT  (M 

The     ICEPLITZ    PRCCEGORE 

OU     150     1=1, N 

Aht  THE  EIGENVALUES  OF  maTkIX  A  =  -D<-C 
IN  [NCKEASING  LnUtk 

01=1/2 

IF  (MODI  I  ,2  )  .i  „.  1  )  AL  =  4DO-200*CCuS(PIN*DI  »  *Ch^ 

IF  (  1  .  I  J. 21  AL«6C0*Ch2 

|»  l  I  ,GT.1.0R.CH2.Gr .DEI  rAlGOTG  120 

A    5PE(   IAi     HKOCEUUKI     FOR       SINGOLAR    CA 

It  |.NI   r.POl  M  )     Gl  Tl    1  12 

POINI     CHAROL 

.     t   ll.JOl 

110    J«l»* 

(AoSiJ-J     I 
Ml,  J)«-.*UO*CJ«S 

1  50 

Ml    nun 

116    J     I  ,  2 
S"0 


114 
116 


lib 


C 

c 

120 


c 

13J 

c 


140 
L5U 

C*«* 
C 

C 


160 


WO 
160 


DO  114  K=1,M 

DKJ=  IAHS(K-J) 

S=S-.5LiO*DK.m ( 1 ,K) 

IFl  J-Eti.  1)  T=S 

I  (  1,  1)=T 

DO  118  J=3,M 

T=F  (  1,  J-ll 

F(  1,J-1)  =  S 

b=2U0*S-F( L, J-2 )-! 

(-(  L,M)=S 

GOTO  150 

END  SINGOLAk  CASE 


COM  INCC 

BE=-500*AL+DSQRT( .2  5D0*AL*AL- l.ODOJ 

BEI=AL-6b 

FUKwAHu    tLlMlNAT  ION 
DO     130     J^2,V 

F( l,J)=F( I, J)+F( I  ,  J-l  )*bEl 
F(  I,f!)=F(l,MJ/(  bfc-oEl  ) 

bALK*AkC  SunSl  1  TOTIJN 
M1=M-1 

DO  140  J=lfMl 
JI=M-J 

F<  I  ,  J I  )={F(  I,  J  I  )♦(-<  I,  Jl  +  1)  l*bEI 
CONTINUt- 

END    GF     TOEPLlTZ    FKGCLOORE. 

INVERSE    FOURIER     IkANSFJkMATI CN    GF    f 

DO  180  I=lfM 
DO  160  J=LfN 
RE( J)=F( J, I  J 

CALL     RFCRT(REtLCG2NtSl,*2, IFERR) 
DO     170    J=lfN 
F( J,  I  J =  KL<  J ) 
COM  INGE 
Kt  TORN 
END 
C 

c 

L  ***********  ***********************  **********************  *****  ***** 

C 
C 

SObKUOTINE  CAPACl 

IMPLICIT  REAL*3(A-H,0-Z) 

LOGICAL  UIR  ,SOL,ELIM,  INFUN 

COMMCN/riLKi/  C<128,128) 

COMMON/BLK^/  F(t^,o1l ,F0N(64,64),RE(64) ,Sl(t>4),Lo{b4) 

CGMMGN/aLK3/  A(12b,J»,X(12d)fw(l^B)tLFx(U8), 

*  CFY1128J  ,LFC(12d) ,HXX( 12 d J  ,HYY{12b>,P(12d),vl(l2d),v"2<l28), 

*  wl(128),SCALES(12tt)tTN(128,2)fIH0(12  6t2)rIPS(128)tID(12  8,4) 
LUMMGN/riLK4/    ACCY,C0NtCH2,0ELTA,H,hSw, IPt  I  0, JO , LGG2N, N, M , 

*  MAX  1  TH ,  I  r  ,DIk,SCL,CL IM, INFON 


C 
C*** 


Trili,  bOJuuOTlNt  CLMPUTfcS  Tht  CAPACITANCE  CATKIX  C  EXPLOITING 
THE  TRANSLATION  IimVAKIANCY  PROPERTY  OF  THE  SOLUTION  OF  HELMHGLTZ. 
EJOATICN  ON  AN  iNFINiTt  STRIP  *ITh  PERIODIC  BOONDARY  CONDITIONS. 


c 


c 
c 


CCJMMMUN    USEU    :     ARRAYS    C ,  «  ,  F  ,  A ,  I  0    4     i  BU  ,     VARIABLES     :    01 R  , DELTA,  N,  I  0  ,  JO 

SC  =  <V.OD0*CH2 

IF ( D  IR  I     GOTO    90 

SINGLE  LAYtk  OPTION, 

00  80  1=1,  I  P 
K=IBO( 1,11 
L= lbDI  1,2) 

HX=HXX(  I  ) 
HY  =  HYY (  I  ) 

KOR=-DSIGN( l.UDO.hX) 
L0R  =  -US1GN(  l.OOCHY) 
Al=  A(  1,1) 
A2=A{  I  ,2) 
A3  =  A( 1,3) 


50 


80 

C 

/  J 

c 


95 


OU  80 
KK=  IhD 
LL=1BD 
K  1  M  A  tJ 
K2  =  I  i-ti 
MK=MIN 
LJ=IA6 
S=SC*F 
IF (0A8 
IAB 
S  -  S  ♦  A  3 
LUNT  IN 
IF (OAd 
b=S+A3 
C( I , J) 
kETUM. 


J=l,  I 
(  J,l  ) 
(J, 2) 

S(K"K 

S(K-K 

0  <  K 1 , 

S(L-L 

1  MK  ,  I 

S(MX) 
SIK  +  K 
*F(Ml 

Ll 

S(HY) 
*f  I  ■•"■. 


K) 

CR-KK) 

N-Kl  )♦  10 

U  f  JO 

J)  *A1*F  (  MN0CK2, 

.Li  .  1.000)  GOTO 

I ^-KK) 

N0(K3,N-K  n  +  10, L J) 


\-K2)*I0,LJ)*A2*HMK,  I AOS ( L-LOR-LL ] ♦ JO ) 

50 


.LE.l.OOO)  GOTO  80 
,  lAbSlL+LUK-LL  l+JOl 


n  N  0  E 

U  n  l  t     L  A  Y  £  H 

.■00     1  =  1,  IP 
>5     1J  =  1  ,  IP 
«( I J)=0.000 

r    I        I  !)()(   1,1) 

Ll  =  Ibi  i  i  ..  ) 

■^■(11 

HY=nYY  I  1  ) 

-     .  I  ;ni  I 
LOR"-OSIGN( 1 

(1,1) 
A  2  "-  m  (  I  ,  2  ) 


liPl  1U\. 


OUt., 
OOCHY  ) 


OU     150    NO"  I  .<♦ 
•     '        1101,102,1      ■  ,'■      ►  1 1 

Ul  •  i  i 

i   I 

I  10 
COM  I 

: 


GOTO  110 

103  CONTINUE 
K  =  KI 

L=LI-LOR 
8  =  A2 
GOTO  110 

104  CONTINUE 
IF(DABSlHX).Lfc.l.ODO)  GOTO  105 
KOR=-KOR 

K=Kl-KQR 
L  =  LI 
6  =  A1 
GOTO  110 

105  CONTINUE 
IF(UABS(HY).LE.1.0D0)  GOTO  160 
LOR=-LOR 

K  =  KI 

L=LI-LOR 

B  =  A2 
110    CONTINUE 
C 

DO  150  J=l, IP 

KK=IBD( J, 1) 

LL=IBD(J,2) 

T=TN(J,1) 

IG=IABS(K-KK) 

S  =  KMINO(  lG,N-IG)  +  IO,IABS(L-LL)+JOI 

KK=ID( J,l ) 

LL=ID( J, 2) 

IG=IABS(K-KK) 

S=S-T*F(MINO<  IG,N-IGH-IO,  I  ABS  (L-LD+JO ) 

KK=IO( J, 3) 

LL=ID( J,4) 

IG=IABS(K-KK) 

S=S*(T-1.000)*F(MINO( IG,  N- IG  )  +  I0,  I  ABS  (  L-LL  )  «-J0  ) 
150    W(J)=W( J)+B*S 
C 
160    CONTINUE 

00  200  IJ=1, IP 
200    C(I,IJ)=W(IJ) 

RETURN 

END 
C 
C 
C ************************************************************** ^^** 

c 

c 

SUBROUTINE  INOUT(NO) 

IMPLICIT  REAL*8(A-H,0-Z) 

LOGICAL  DIR, SOL, ELIM, INFUN 

C0MM0N/BLK2/  F ( 64,64 ) , FUN( 64, 64) , RE ( 64 ) , SI ( 64 ) , LG ( 64) 

COMMON /BLK3/  A(12  8,3),X(128),w(128) ,CFX<128), 

*  CFY(128),CFC(128),HXX(128),HYY(128),P(128),V1(128),V2(128), 

*  Hi  (  128 ),SC ALES ( 1 28 ), TN( 128, 2), I  BO (128,2), I  PS (128), ID (128, 4 1 
COMMON/ 8LK4/  ACCY ,C0N,CH2 , DELTA, H, HSQ , I P , 1 0 , JO, L0G2N, N, M , 

*  MAXITR, IT, DIR, SOL, ELIM, INFUN 

COMMON/F/  FMTK  10) ,  FMT2  (  10  ) ,  FMT3(  10>,FMT4<10) 

OAT  A  ATEX T,BT EXT, CTE XT, DTE XT /'NEUMANN', 'DIRICHL.' ,'DU/DN=G,f 'U=G'/ 

IF(DIR)ATEXT=BTEXT 

■ 


IN  THIS  SUBROUTINE  ALL  THE  USER'S  SUPPLIED  INPUT  IS  COLLECTEO. 
THE  INPUT  DATA  IS  PRINTED  OUT  AS  A  CHECK. 

GOTO  (  100, 200, 300, 400, 500), NO 
CONTINUE 
A  NUMbER  OF  PARAMETERS  MUST  BE  INITIALIZED. 
THEY  ARE: 

CON  -  THE  CONSTANT  IN  THE  HELMHOLTZ  EQUATION. 
LUG2N  -  AN  INTEGER  SUCH  THAT  N=2**LOG2N  IS  THE  NUMBER  OF  MESH  POINTS 

ACCROSS  THE  STRIP 
IP  -  THE  NUMBER  OF  IRREGULAR  MESH  POINTS. 
M  -  THE  NUMBER  OF  MESH  POINTS  USED  ALONG  THE  STRIP.  ITS  DEFAULT 

VALUE  IS  N. 
IH  -  AN  INTEGER  SUCH  THAT  2*IH  IS  THE  NUMBER  OF  MESH  SIZES  ACROSS 

THE  REGION.  IH  MUST  BE  LESS  THAN  N/2.  ITS  DEFAULT  VALUE  IS  3N/8. 
DIR-  EQUAL  TO  .TRUE.  FOR  THE  DIRICHLET  PROBMLEM, 

EQUAL  TO  .FALSE.  FOR  THE  NEUMANN  PROBLEM 
SOL  -  EQUAL  TO  .TRUE.  IF  THE  ANALYTICAL  SOLUTION  IS  KNOWN, 

.FALSE.  OTHERWISE. 
ELIM  -  EQUAL  TO  .TRUE.  FOR  THE  GAUSSIAN  ELIMINATION  OPTION, 

EQUAL  TO  .FALSE.  FOR  THE  CONJUGATE  GRADIENT  OPTION. 
INFUN  -  EQUAL  TO  .TRUE.  IF  AT  EACH  CALL  OF  THE  SUBROUTINE  DEKEN 
NEW  DATA  F  IS  READ  IN  ,  .FALSE.  OTHERWISE. 
THE  DEFAULT  VALUE  OF  INFUN  IS  .FALSE. 
MAXITR  -   THE  MAXIMUM  NUMBER  OF  RIGHT  HAND  SIOES. 

ACCY-THE  REQUESTED  ACCURACY  OF  THE  CONJUGATE  GRADIENT  SUBROUTINE. 
DELTA  -  THE  MACHINE  PRECISION 

IF  POSSIBLE  THE  DOMAIN  SHGULD  BE  ORIENTEO  SO  THAT  M>=N. 

10, JO  ARE  THE  COORDINATES  OF  THE  POINT  CHARGE  USED  IN  SOLVE! .TRUE. ) . 

THEIR  VALUES  SHOULD  NOT  BE  CHANGEO. 

READ(5,1)  FMT1,FMT2,FMT3,FMT4 
1       F0RMAT(10A8) 

KEAD(5,FMT1)DIR,ELIM, INFUN , SOL, ACCY , CON, DELTA , LOG 2N , MAX  I TR, IP 

WRITE(6,FMT1  ) DIR, ELIM,  I NFUN, SOL , ACC Y , CON, DELT A , L0G2N, MAX  I TR , I P 

IH=3*2**(L0G2N-3I 

N=2**L0G2N 

M  =  N 

10=2 

J0  =  2 

H=1.000/DFL0ATl IH) 

HSQ=H*H 

CH2=C0N*HSQ 

RETURN 
C 

200    CONTINUE 

C***     INPUT  THE  REUUIRED  INFORMATION  ABOUT  BOUNDARY  OF  THE  REGION. 
C***     IBJ  CONTAINS  THE  COORDINATES  OF  THE  IRREGULAR  MESH  POINTS. 
L***    DISTANCES  (IN  MESH  SIZES)  FROM  THE  IRREGULAR  MESH  POINTS 
C***    TO  THE  BOUNDARY  CURVE  IN  THE  I  A  J  DIRECTIONS  ARE  STORED 
C*«*     IN  HXX  A  HYY,  RESPECTIVELY. 

C**»    TN  CONTAINS  THE  TANGENTS  LF  THE  NORMAL  DIRECTIONS. 
C 

C***    THE  IRREGULAR  MESH  POINTS  CAN  BE  ORDERED  IN  AN  ARBITRARY  WAY. 
C»*«    HOR  EACH  IRREGULAR  MESH  POINT  THE  USEK  MUST  SUPPLY  THE  FOLLOWING: 
C+**  ItJtHl.HJ.TG   FUR  THE  DIRICHLET  PROBLEM, 

L***     I  , J, HI ,HJ, TA.TB   FOR  THE  NEUMANN  PROBLtM. 
C***  COURDINAIES  I  A  J  (  BETWEEN  2  AND  N-l  FOR  I,  BETWEEN  2  AND  M-L  FOR  J) 


I-DIRECTION  IS  ACROSS  AND  J-DIRECTION  ALONG  THE  VERTICAL  STRIP. 

NOTE  THAT  THE  VALUES  OF  1  ANO  N  FOR  I  AND  1  AND  M  FOR  J  ARE  NOT  ALLOWED. 

NO  POINT  (I, J)  WITH  THREE  NEIGHBOURS  ON  THE  BOUNDARY  IS  ALLOWED. 

HI  A  HJ  -  SIGNED  DISTANCES  IN  MESH  SIZES  FROM  THE  POINT  (I, J) 

TO  THE  BOUNDARY  CURVE  ACROSS  ANO  ALONG  THE  STR IP ,R ESPECT IVELY. 

THE  ABSOLUTE  VALUES  OF  HI  AND  HJ  MUST  BE  BIGGER  THAN  DELTA. 

IF  A  POINT  (I, J)  HAS  ONLY  ONE  NEIGHBOUR  ON  THE  BOUNDARY  CURVE 

SET  THE  OTHER  VALUE  OF  HI  OR  HJ  TO  ANY  VALUE  LARGER  THAN  ONE 

IN  ABSOLUTE  VALUEIwITH  THE  CORRECT  SIGN). 

THE  VALUE  OF  ONE  ADDITIONAL  VARIABLE  TG  IS  NEEDED  IN  THE  DIRICHLET  CASE 

AND  THOSE  OF  TWO  VARIABLES  TA  AND  TB  IN  THE  NEUMANN  CASE. 

TG  -  THE  TANGENT  OF  THE  SMALLER  OF  THE  ANGLES  BETWEEN  THE  NORMAL 

THROUGH  THE  POINT  (I, J)  AND  THE  AXES  ACROSS  AND  ALONG  THE  STRIP. 

TA  -  THE  TANGENT  OF  THE  ANGLE  BETWEEN  THE  NORMAL  AT  (I-HI,J) 

AND  THE  AXIS  ACROSS  THE  STRIP. 

TB  -  THE  TANGENT  OF  THE  ANGLE  BETWEEN  THE  NORMAL  AT  U,J-HJ) 

AND  THE  AXIS  ALONG  THE  STRIP. 

THE  VALUES  OF  TA,TB  AND  TG  SHOULD  ALWAYS  BE  NON-NEGATIVE. 

IF  A  POINT  (I, J)  HAS  ONLY  ONE  NEIGHBOUR  ON  THE  BOUNDARY  CURVE 

SET  THE  OTHER  VALUE  OF  TA  OR  TB  TO  ANY  VALUE. 

WARNING  : 

THERE  IS  NO  TEST  FOR  POSSIBLE  BUGS  IN  THE  BOUNDARY  REPRESENTATION. 

READ  IN   I,J,HI,HJ  AND  TG  FOR  THE  DIRICHLET  PROBLEM 
OR  TA,TB  FOR  THE  NEUMANN  PROBLEM. 

IF(DIR)  GOTO  211 

DO  210  1  =  1,  IP 
READ(5,FMT2)  IBD( 1 , 1 ) ,  I  BD( I , 2) ,HXX(  I ) ,HYY( I ) , TN( I , 1 ) ,TN( I,2> 
wRITE(6,FMT2)  IBD(  1 ,  1 )  ,  IBD(  1 ,2  )  ,HXX(  I )  ,HYY< I ) , TN( I , 1 ) , TN U ,2 ) 
CONTINUE 

GOTO  215 

CONTINUE 

DO  212  1=1, IP 
READ(5,FMT2)  IBD( 1 , 1 ) ,  IBD( I,2),HXX( I),HYY< I), TNI  1,1) 
WRITE(b,FMT2)  I80( 1 , 1 ) ,  IBD(  I  ,2 ) ,HXX(I),HYY( I),  TNI  1,1) 

CONTINUE 

CONTINUE 

RETURN 

CONTINUE 
THE  USER  SHOULD  SUPPLY   N*M  VALUES  OF  THE  FUNCTION  F(X,Y)  AT  EACH 
MESH  POINT  OF  THE  IMBEDDING  RECTANGLE. THE  VALUES  AT  THE  EXTERIOR 
MESH  POINTS  CAN  BE  CHOSEN  ARBITRARILY. 

IF(.NOT.INFUN.AND.IT.NE.i)  GOTO  320 
REAO(5,FMT3)  I (FUN( I , J ) ,  1  =  1, N) , J=l, M ) 
WRITE(6,FMT3)  <(FUN(I,J),I=1,N),J=1,M) 
320    CONTINUE 
RETURN 

CONTINUE 
FOR  EACH  SET  OF  DIFFERENT  DATA,  I.E.  MAXITR  TIMES, 
THE  USER  SHOULD  PREPARE  THE  FOLLOWING  DATA  : 
FOR  EACH  IRREGULAR  MESH  POINT  VALUES  OF  VIII)  ANO  V2(I). 

VI  (V2)   IS  THE  NEUMANN  OR  DIRICHLET  DATA  AT  THE  POINT  ON  THE  BOUNDARY 
CLOSEST  TO  THE  IRREGULAR  MESH  POINT,  IN  THE  DIRECTION  ACROSS ( ALONG) 
THE  STRIP.  IF  IN  ONE  DIRECTION  THE  BOUNDARY  IS  AT  A  DISTANCE  LARGER 
THAN  H  AN  ARBITRARY  VALUE  CAN  BE  ENTERED  FOR  THE  CORRESPONDING  VI  OR  V2 
THE  BOUNDARY  MESH  POINTS  ARE  ORDERED  AS  IN  THE  SUBROUTINE  BOUNDR. 

DO  410  1=1, IP 

88 


READ(5,FMT4)  V1(I),V2UI 
WRITE(6,FMT4)  Vl(I)tV2(I) 

410    CONTINUE 

RETURN 
C 

500    CONTINUE 

C***    IF  THE  LOGICAL  SOL=.TRUE.  THE  VALUES  OF  THE  NEGATIVE  LOGARITHM 
C***    OF  THE  ABSOLUTE  VALUE  OF  THE  ERRCR  ARE  PRINTED. 
C 

HRITE(6,6  50)ATEXT,CTEXT 

WRITE(6,620)C0N 

WRITE(6,630)N,M 

*RITE<6,640) IP 

IF  (SOL)  GOTO  545 

WRITEI6,542) 
542    FORMATC  SOLUTION*  I 
C 

C**    AOO  PRINTOUT  IN  REQUIRED  FORMAT 
C 

GOTO  555 
545    MRITE(6,660I 

IU=0 

T  =  000 

V=ODO 

DO  550  J1=1,M 

J=M*1-J1 
C***    IF  SOL=.TRUE.  THE  USER  SHOULD  PREPARE  N*M  VALUES  OF  THE  SOLUTION 
C***    U(X,Y)  AT  EACH  POINT  OF  THE  IMBEDDING  RECTANGLE. 
C 

REAO(5,FMT3)  ( RE ( I ) , 1=1 , N ) 

DO  548  1=1, N 

S  =  F(  I  ,J)-RE( I ) 

S=DABS(S) 

LG(I )=DMAX1(OCO,-DLOG10(S*1D-20) )*.5D0 

IF(S.GE.1D2*ACCY)  GOTO  548 

V=UMAX1( VtS) 

T=T*S*S 

IU=IQ*1 
548    CONTINUE 

WRITE(6,700)(LG( I  ),I  =  1,N) 
550    CUNTINUE 

IF( IQ.Nt.O)T=DSQRT(T/DFLOAT( 10) ) 

MRITE(6f 11001 

WRITE(6, 12001  IQ,V,T 
555    CONTINUt 

RETURN 
620    FORMATC  C=',1P010.3) 

630    FORMAT! •  NUMBER  OF  POINTS  IN  RECT .R EG  ION= • , I  3 , • ♦ • ,  I  3 ) 
640    FORMATC  NUMBER  OF  BOUNDARY  P0INTS=',I3) 

650    FORMATC  Iht  IMHOLTZ  tQN.NlTH  ',A8,'  BOUNDARY  CONDITIONS'/ 
*•  LU*CU=F  I*  REGION  •/ 
••  ',A8,'  ON  BOUNDARY'/  ) 
660    FORMAK/'  -L OG < I E RROR | ) ' / ) 
700     KJRHAT ( 1X.64I2) 

L100   FORMAT!//'  NO.  POINTS  MAX  A  L2  -  NORMS  OF  THfc  ERROR') 
1200   FORMATC 15, 8( 1P015.5) ) 

END 
C 
C 
(^•♦••♦•♦♦•♦♦♦♦♦•♦♦•♦♦♦••♦♦♦•♦♦♦♦♦•♦♦♦••♦•••♦♦♦♦♦♦•♦•♦♦••♦♦•*******« 

C 


c 

c*** 
c*** 
c 


10 


20 


30 


bO 


60 


70 


SUbRUUl INE     S>MMLG(N,X  ,  ii,P,  VI  ,  V2 ,  »  ,  EPS  MAC  ,  ACC  Y  .  I  NT  MAX,  I  S  TOP  ) 

IMPLICIT    RFAl*ti     (A-h,C-Z) 

DIMENSION    X  (N)  ,  i  I :. )  ,  P(N)  ,  VI  IN),  V2(N),w(N) 


CuNJGOATE    GRADIENT    SUBROUTINE    uY 
STANFORD    RFPUKT    C S- 399 , NOV . 1 97 3. 


C.C.PAIGE  AND  M. A  .S AUNOLR S , 


EPS=ODO*EPSMAC 

CALL  AT 1MESIX.P ,N) 

UU  10  1=1, N 

Vl(  1  )  =  B(  1)-P<  I) 

CUNT INCE 

CALL  DNORMI VI ,\,EPS,D1) 

UKN0RM=D1 

CALL  ATIMESI VI ,P,N) 

ALPHA=ODO 

DU    20     1=1, N 

hi  I  )=V1<  I) 

ALPHA=  l/KI  ).*P(  I  J  +  ALPnA 

CONTINUE 

Du     30     1  =  1, N 

V2(  I  )  =  Pl  l)-ALPHA*Vl I  1 ) 

CUNTINOE 

CaLL    GNCkM(V2,N, EPS, BETA) 

GbAK=ALPHA 

DBAR=BETA 

D2=0D0 

DNnRMX=OUO 

GiMU^A=UABS(ALPt-A)+BETA 

EPSA=DNGKMA*EPS 

EPSX=EPSA 

ITN  =  0 

ISTCP=0 

DLwNOR=GSWRTt Ll**2+D2**2) 

CGNGRM=3RNGRM*BETA/  (OABblGBARJ+EPSA) 

BESTNM  =  DMIM  ( DLGNCR tCGNGRM ) 

IK  ITN.EC.  ITNNAX)     ISTCP=3 

If-(BESTNM.LE.EPSX)     ISTGP=2 

IF ( otSTNW.Lt.ACCY  )  ISTCP=1 

IF(  1  STOP.NE.O)  GOTO  100 

*KlTE(6,  1010)  lTNtX( l)t OLGNORfCGNORM 

CALL  ATIMES(V2,P,N) 

ALPHA=COO 

DO  60  1=1, N 

ALPMA  =  v2(  I ) *P(  1  J  +  ALPHA 

CGNTINGE 

00  70  1=1, N 

T  =  V2(  I  ) 

V2(  l)  =  H(l)-ALPHA*1-BtTA*Vl(  I) 

Vl(  1)  =  T 

CONTINUE 

GL0B=6ETA 

CALL  DNGRMI V2 , N , E P S A , BE T A  ) 

GAMMA=GSQRT ( GB AP **2 +0L0 B** 2 ) 

L5=GBAR/GAMMA 

SN=OLDB/GAMfA 

DELTA=CS*DBAR+SN* ALPHA 

GBAR=SK*CBAK-CS*ALPhA 

EPSLN=SN*HETA 

DBAR=-CS*BETA 

URNGkM=SN*gRNGKM 

Z  =  01/GA^A 


90 


do 


100 


11J 

120 


101G 
1020 


S=Z*CS 

T=Z*SN 

00  80  1  =  1, N 

X(  I )=<  w(  I  )*S+V1  (  II*D+X(  I) 

Hi   I  »  =  W< I )*SN-V1 ( 1  )*CS 

CUNT INOE 

S=UL OB *DABS( ALPHA)* BETA 

IKONORMA.LT. S)  ONORMA=S 

ONGRMX  =  Z**2*0NG KM  X 

EPSA=UNCRMA*CPS 

EPSX  =  DS(jKT(  CNGkMX)*E?SA 

D1=G2-GELTA*Z 

02=-EPSLN*Z 

[TNMTN+1 

GUTC    5  0 

it-  (OLUNOR.LE.CGNORM)  IST0P=-ISTUP 

ik  istcp.lt  .0)   goto  120 
Zba;<=oi/gbar 

do     110     1=1, N 

X 1  1  )  =  m (  1  )*ZbAK+X(  1  ) 

CUNT INUE 

WKITLIO.IOIO)     ITN.X(l)  ,  OL(.NCR,CGNGRV 

1020)     ITN,  ISTGPtACCY,EPSX,BESTNM 
CGNGkfrtDLQNURtCiRNORM 


wRl  TE(6 
MKITE(6i 1030) 
RETURN 
FORMAT! 


I  d,  1PD2C.  10, 1P6015. 5) 


FORMAT (  •  1*  , 


,  16X.U0, 


,  1BX,  I  10, 


Uiu 


*  /  ,  •  Nl.  Of   ITT  h ATIGNS: 

*  /,'  STUt>PlNu  CONDITION  WAS 

*  /,'  NURM  UF  RtSIGJAL  wAS  -UJOIRED  TU  Bt:  'vlP015.5f 

*  /,'  ESTIMATE  LF  KEASUNABLE  \CkM:  »tiP0l5.5i 

*  /,'  I:>1I*AIF  Of  NORM  ACTUALLY  OBTAINED:   ',1^015.5) 
FORMAT  I  «0E ST  1KATES  CF  NGRM  gf  final  residual:', 

*/,'  L J:' .1P015.5/ , '  LJ  INCOMPLETE :' ,1PD15. 5/, •  QR: • , 1PD15.5) 
fcND 


C 

C 

c**«* 

c 

c 


***  *  *  *  *  *********  *  *  «,M«M,M«*  *******  ******  *******  ************ 


10 


20 


C 

c 
(_•♦  ♦♦ 

L 

c 


SUB 

Imp 
DIM 
S  =  U 
OU 
S-  V 
GUN 

If  ( 

OU 
*(  1 

GUN 
R|  1 
ENO 


L  IC 

ENS 
CO 

10 

( I ) 

T  IN 

A  =  lj 

OET 

GO/ 

20 

»»V 

T  IN 


I  I.I  CNLRMV.N, LPS, BETA  ) 
IT  RFAL*8  (A-(-,l!-Z) 
ICN  V(N) 


1  =  1  ,N 
**2*S 
UE 

(  S) 

a. i  r  .lps  ) 
BETA 

1  =  1, N 
(  I  )*S 

OF 


*LIa=lp:>*.500 


*****«**««*****i«<.******»*******»******t**»****»************** 


iUtikUUl  INE     Al  1  v,LSU,PP,NN) 
IMPLICIT    REAL*MA-H,U-Z  ) 
UI MENS  ION    V (NN)  fPf  INN) 
CUMMON/bLKl/    (.(128,128) 
CJMMON/AT/     Pl(12b) 


C 

L*** 
C*** 
C 

c*** 
c 


u 

20 

c 

c*** 
c 


30 

<*0 


c 

c 

c***< 

c 

c 


c 

L*** 
C*** 

c*** 
c 


CUMPUIES  THE 
COMPGN  USEO 


PRODUCT 
AKPAYS 


UF  MATRIX 

C  ANU  PI. 


C'C  AND  VECTOR  V. 


P  1  =  C  *  V 


NN 


DU  20  1=1 

S=u.000 

DU  10  J-=ltNN 

S=S*C( I  ,  J)*v< J) 

Pl(  I  )  =  s 

P=C  'CV 

DU  40  J=1,NN 

S=0.000 

DU  30  1=1, NN 

S=S+C(  I  ,  J)*P1 ( i  ) 

PP(  J  )=S 

RETURN 

END 


$$*^<-^<.  **************:;:******<•********************************* 


SUBROUTINE  CCCU^PI  SCALES, IPS,  IP) 
IMPLICIT  HtAL*6  IA-htO-Z) 
DIMENSION  SCALES!  IP)  ,  IPS(  IP) 
CUMMUN/dLKI/  C(12b,12d) 


1 
2 

5 


ID 


LU  DEC0MP0S1SICN  UF  MATRIX  C 
THE  RESULT  LU  OVERWRITES  THE 
COMMON  USFD  :  MATRIX  C. 

DU  5  1  =  1,  IP 

IPS!  I  )  =  l 

RQwNRM=0D0 

DU  2  J=l,  IP 

IM  f<CfcNRM-DAbSU.(  I,  J)  1)      l»2»2 

KUWNKM=DAbS(Cl I , J) ) 

CONT  INUE 

lF(KUwNRW)     J, A, 3 

SOALLSl  I  )=1D0/KLwNRM 

OuTu    b 

«KlTt (b, 100) 

SCALESl I )=OD0 

CUNT INOF 

IRM1=1P-1 

DU     1/    K=1,IPM] 

blO=UDC 

DU     11     l=K,IP 

IS=  IPS(  I  ) 

SiZE  =  DAdS(C(  1S,K)  )*SCalF.S(  IS) 

IMSlZfc-iilG)     11,11,10 

blG=SIZE 

IDXPIV=I 


At-TEK    FURSYThfc-MULER. 
INPUT    MATkIX    C. 


92 


11 

12 

13 
14 

15 


16 
17 


1  J 
19 
I  JO 

C 
C 

c 
c 


c 

t  *** 

c*** 
c 


CUNT  INUE 

IF  (BIG)  13,12,13 

wRITE(6,100) 

OUT  [J  17 

IF (  IDXPIV-K )  14,  lb,  14 

JMPS(K) 

IPS(K) =  IPS< IDXPIV) 

IPS(  IDXPIV) =J 

i\P  =  I PS(K) 

PI VUT=C( KPtKI 

KP1  =  K*  1 

DO  16  1=KP1,IP 

IS=IPS( 1 ) 

LM  =  -C(  IS,K) /PIVCT 

C( IS,K)=-EM 

DU  16  J=KP1,IP 

C(  IS,J)=C(  I S, J )+tM*ClKP,  J) 

CONTINUE 

CUNT INLE 

*P  =  I  PS(  IP) 

IK  (C(KP,  IP)  )   19,  Id,  19 

MRtTE(6f 100 ) 

RETURN 

FukMATJ/'  SINGULAk  MATRIX1) 

END 


<t*<-<:<;^*«:«:<:*<:*<:<:t<=<^*^^'?**t***<!t*<-***^*«:«*<!****«*****<r********<r 


3 
4 


SUbhOUT INE     fVLCSIX,*, SCALES, IPS, IP) 
IMPLICIT     PLAL*H     (A-H,u-Z) 

DIMENSION    X(IP),*UP), SCALES!  IP),1PS(1P) 
ClJMMCN/ttLKl/    (.(  l2o,128) 

SULV1NG   Lu  DECOMPOSED  LINEAR  SYSTEM  UF  FuNS. 
THE  SOLUTIUN  IN  a  ,  TnE  RHS  IN  k. 
CO^MLN  USED  :  NATKlA  C. 

IP1=IP*1 
!:>=  IPS(  1  ) 
XI  I )=*( is) 
DU    2    1=2,  IP 
is=  ips( i ) 

IMU  I-  1 

SUM^ODC 

JU     I     J  =  I  ,  I  M  I 

SUM«SUP*CI I \, J ) * x ( J ) 

A(  I  l*tol  IS)-Sum 

IS=IPSI II I 

a l IPI-XUPI/Cf I  St  IP] 

DO    4    ! mACk=2, IP 

|*IPl- IHACK 

IS=  IPS!  1  ) 
IPL 1=1 »J 

sum^odo 

UU     3     J"  IPL  I  ,  IP 

SUM=  Su"*U  I  S, J  I  «X ( J  ) 

■  i  I  )  -  I  X|  1  ) -SUPI /C(  I  S,  I  ) 

- t    TUKN 
END 


C 

c 

c 
c 


c 
c**» 

L*** 

c 

1 


5 
6 
10 


20 


30 


40 

c**« 

60 


«,*<<«<  *<.<■«<'**<■***  t*****************1**************************** 


SUom.UT  INL     kFI:kHA,M,S,IFS»IFERK) 
IMPLICIT    RfcAL*tt(A-hfO-Z) 
DIMENSION    A( L  )  ,S(  I) 

CGuLtY'S    REAL     J-CUKltK     TRANSFORM    PRGGkAM 

USING    FFT     SUBRCUTINF     -    FORT 

rflTH    ChANGFS     ALLOWING     THE     STGkAbt    <JF    ALL     FUUkltR    COEFFICIENTS 

IN    A    VECTOR     (l:Nl     :     SIN    <i    CUS    AMPLITUDES    CF    GROWING    FREQUENCY 

1=1    Tl    N/2-1    AKb     Iin    PUS.2*I+lf  2*I+2;COS(OJ     IN    PuS . 1 , COS ( N/2 )     IN    PGS.2 

lFEkkS=0 

N=2**M 

NV2=N/2 

NV4M1= N/4-1 

MMl=M-l 

IF  (  l«bS(IFS)-l  )  10,lO,t> 

IF ( MP-N)b,20,20 

[f  I  KKS=1 

NP  =  N 

MP=M 

CALL  FGkT( A,M, StO, IFEkRl  ) 
[ F I  RRS=IFERRS+IFFRR1 
K  D=  N  F  /  N 

kt  =  ko 

i\PV4  =  NP/4 

IF  (  IFS)30,b0,6C 

CALL  FLKTI A  ,MMl,St-2, IFERR2) 

lFLkRS=  IFERRS+IFLHR2 

<*G  K  =  1,NV4M 
J=NV2-K 

A1R=A<  2*K  +  l  )  +A(  2*J+1) 
Al I=A( l*K+l )-A(k*J+2J 
A2K=A( 2*K  +  2 )+A(2*J+2) 
A2I=A(2*J+1 )-A(2*K+l) 
KKT=NPV4-KT 

AHR=A2R*S(KKT)+A2I*S(KT) 
AwI=A2l*:>(KKT)-A2k*b(KT) 
A(2*K+1)  =  (  A1R  +  AWR)/4.01)J 
A(  2*K+2)  =  (A11+Ak»n/4.0U0 
A(2*J+l)=<Alk-AWRI/4.0D0 
A(ii*J+2)  =  (AWl-All  J/4.OD0 
KT=KT+KO 
T=A( 1 ) 
A(l)=( T+A<2) ) /2.UU0 

ORIGINALLY     A(2)  =  A(N«-2)=0    AND 
A(2        )=<T-A(2) 1/2.000 
A(NV2*l)=.5O0*A (NV2+1) 
A(NV2*2)=-.5U0*A(NV2*2) 
lFEkk=  IFFRkS 
RETURN 

DU     60    K=lfNV4Ml 
J=NV2-K 

A1R  =  A( 2*K+1  )  +  A(2*J*-  1  J 
Al l=A< 2*K+2)-A(2*J*2J 
A*k=A(  2*K*l  )-A(  2*J+1) 
A*  I=A(  2*K*2)«-A<  2*J+2) 


A(N*1)=(T-A< 2) )/2 


94 


8<J 


KKT=NPV4-KT 

A2K=A*k*S(KKT)-AwI*S(KT) 
A2l=AWk*S(KD+AWi*S(KKTJ 
AU*K+1)=AIR-A2  I 
A(2*K«-2)=A1 i*A2k 
A<2*J+  l)=AiK*A2 I 
At  2*J«-2)=A2R-A1 I 
KT=KT+KD 
T  =  A(  1) 

UKIG1NALLY  £=A(N+1) 
Z=A(2) 
A(  1  )=T*Z 
A(2)=T-Z 

A( NV2+l)=2.OD0*A(NV2+U 
A(NV2+2)=-2.OD0*A(NV2+2J 
CALL     hL'kTIA,  M*l  ,S,2  ,  IFtkk2) 
It-LKkS=IFEkRS*IFDHK2 
GUTG    50 
LNU 
C 

c 

C 

c 

SUBHUUTINE    F0kT(Af y.S.IrSt IFfcRK) 

IMPLICIT    REAL*8(A-Hf0-ZJ 

lJ  I  M  E  r>i  S  1 0  N    A(  I  )  ,  S(  1  ) 

[F(M)£t2v5 

IFEPR=1 

RETURN 

[FERR=0 

N2=2*N 

It-  (  [A6S(  IF  S  )  -  i  )  ^00,2uO,  10 

It-  (  N-NP)20f  20il2 

It-  E  rK  =  l 

buTu    ^CO 

NV2*N72 

NM1  *N-  1 

J»l 

Du     jO     l«ltNMl 

IF ( I  .GE.JJGC Tu    22 

A(  t* l-l )=At  2*J-  1  ) 
A  ( 2  *  J  - 1 )  «  T 
r*A(2*l] 

A12*  I  )=A(<?«J) 
A(  2*J)  =T 

K  =  NV<> 

It  ( K.GE  .J)GOTG    30 

J=  J-K 

K  ■  K  /  2 

GOTO    2<. 

J- j  •  - 

If  (  IFSI       .     »3< 

FN»N 

UU      3<.      I  =  1  .  N 

I  |   ■ -  I  -  I )    M 2*1-1) /FN 

A(2*U— A(2*1)/FN 

ou  t  o   i«  1 1  n  , ; 

T  -  A  (  2  ♦  I  -  1  J 


2 
1 
5 


10 
12 

20 


2o 


11 
2<» 


10 

3<! 


46 


Ai2*l-ll=T+A(2*l+lJ 

A(2*I+1 )  =  T-A< 2*1+1) 

T=Al 2*  1  ) 

A(2*l  )=T+A<  2*1+2) 
tu     A(2*I+2)=T-A(2*I+21 

IF (M-l  )2t It  50 
50     LEXP1=2 

LEXP=8 

NPL=2**MT 
oO     DO  130  L^2,P 

0U  bO  l=2,N2,LLXP 

U  =  I+LEXP1 

12  =  1  1+LEXP1 

!3=I2+LtXPl 

T  =  A(  1-1  ) 

A(  I-1)=T+A(  12-1  ) 

A(  12-1  )=T-A(  12-1) 

T=Al  I  ) 

A(1)=T+A(12) 

Al  12  )  =  T-A(  12) 

T  =  -A( 13) 

T  l  =  A(  13-1  > 

A(  13-1  )=A(  1  l-l  )-I 

A(  I3)  =  A(  II  )-TI 

A(  I  1-1  )=A(  1  1-1  )  +  l 
60     A(  I1)=A(  11)  +TI 

IFIL-2)  120,  12C, SO 
90     KLAST=K2-lfcXP 

Jj=iSiPL 

00  11C    J=4,LEXPl»2 
NPJJ=NT-JJ 
Ok=SINPJJ) 
UI=SIJJ) 
1LAST=J+KLAST 

Ou     100     I=J t  ILAST,LEXP 

1  1=  1  +LEXP1 
12=1  1+LEXP1 
U  =  I2+LEXP1 

T  =  Al  12-1)*UR-A(  I2)*Ul 

TI=A( I2-1)*U1+A( 12)*UR 

Al  12-1  )=A{  1-D-T 

A  I  I  2  )  =  A  <  I  )  -  1  I 

A(  1-1  )=A( 1-11+T 

Al  I  >=Al  I  )+T  I 

T=-A(  I3-1)*U1-A( I3)*0R 

T  I  =  A( I  3-1) *UK-A(I 3)  *UI 

Al  13-1  )=A(  I  1-1  )-T 

Al  13)  =  Al  U  >  — T I 

Al  1  1-1  )=A(  1  l-l  )+T 

100  Al  I  1  )  =  A(  II  )  +U 

110  JJ=JJ+NPL 

120  LEXP1=2*LEXP1 

LEXP=2*LtXP 

130  NPL=NPL/2 

L40  IF  I  IPS  )  145,2,1 

145  00     150     1  =  1, N 

150  A12*I  ) =-A(2*I  ) 

160  GuTU     1 

POO  NP=N 

MP  =  M 
NT=N/4 


96 


205 
210 


220 


2jO 


MT=M-2 

IF(MT)260,260,2C5 

THE1A=.7  85  3  981633974483D0 

JSTEP=NT 

JDIF«NT/2 

S( JC1F  J  =  LSIM ThETA) 

If-  (MT-2)  260,220,220 

DU    250    1  =  2, l"T 

ThL  TA  =  TH(rTA/2.0CO 

JbTCP2=JSTEP 

JSTEP=JDIF 

JUIP=JDlF/2 

SUDIf  J=LSlMThETA) 

Jti=NT-JOIF 

S(  JC1)  =DC0S(  THt  TA) 

JLAST=NT-JSTEP2 

IF(JLAST-JSTEP)*:5C,2  30,230 

DU    240       J=JSTEP,JLAST ,JSTEP 

JC=NT-J 

JU= J+JLlf 

S(JO)=S(J)*S(JCl)*b(J0IF)*S(JC) 

CUNT  INUL 

IF  (  IFSJ20, ^61,20 

re  fuk.\ 

END 


Z<*0 
250 
260 
2t>l 

C 
C 

c 
c 

INPUT     OIKIChLtT     DATA    FOR     THE    REGlUN     IN    FIG. 2. 
N*i6f     IP=27,     U=lX*X*-Y*Y-ll/4»    F«-l. 

FuKMAT     SIAUMNTS 

|4L2t3Dl4.5i 315) 
12I5.4F15.10J 

<<»F  15.  10) 
(2F  15. 1U) 


I   .1  I  i    .i 
T     F     F     T 

I 


(mkA^ETERS 

1  .oo-t 

HI 


0.000 
hj 


I. 00-15 
TG 


1  27 


13  6  -  l.U 

I  i  •/  -0.91  bu/w  run 

II  10  -0.656d5<«/^5 
l)  11  -0.1961524227 
12  l .  -0.4721359 

li  I  :)  -0.3166247 

10  11  -I.  H6b<:W904 

-  l  >  -2.3 16624  7504 

8  1  j  1.3       » 3  5  5  C  9  5 

1  it  .302435!      • 

/  12  0.7562CV 

6  LI  .... 

b  10  .        '   !    •  ^  1  <.0  S  1 

6  .32323909<H 

u  b 

b  /  .32  3239 


-1.7320508076 
-2.31662479  14 
-1  .  316b/". 

-  ).  3166247904 
-0  .  4  72 1 3  595  50 

-  J.  196152^,22  7 
-0.65685424 • 
-0.916079/831 
-1.0 

-0 • 43281 18 1 77 
-1.4328118177 
•0.11  .203 

-1.172     I     i 
-2.172U 
(.172     •     • 
2.1720304 


u.o 

0. 2003000000 

O.^UuCOOOOOO 

3.6000000 

l.OuJUJJOOO 

0.600COOJOOO 

.  .  )00000000 
0.20000  I 

. 
0.771 

.    .168502  rSl 

0.231  M 

.  1 

0.07  7lOb28<.<. 

. 
J.  0/  H 


6 

6 

0.22  144146S1 

1 

.  1  7203C4203 

0. I542l256dd 

6 

5 

C. 0405242848 

0 

.17203C4203 

0.2313188532 

7 

4 

0.7562036828 

1 

.43281  18177 

0.6168502751 

I 

1 

0.  3024355uS5 

0 

.432511817/ 

0. 7710628438 

8 

J 

-1.  732050bo/6 

1 

.0 

0.0 

9 

J 

-  0  .  7  3  2  0  5  0  B  o  7  o 

0 

.  732G5C8076 

1.0000000000 

9 

* 

-1.0 

1 

.  7320*>Oo076 

o.O 

10 

5 

-0.2679491924 

1 

.0 

0.5OOC0OOOOO 

11 

o 

-1.0 

J 

.26  /9491924 

0. 50U0OO0OO0 

12 

7 

-1.73205Ofc076 

1 

.0 

0.0 

13 

7 

-J. 7320506076 

0 

.7320503C76 

l.OOOCOOOOOO 

V  11  1  ) 

-0.0000000000 
-0. JOOOOOOOOO 
-u. J  100000000 

-o.ooooococoo 

-o.oooooooooo 

-O.OOOOGCOCCC 

-o.oooooooooo 
-o.ooouooocoo 

-u.  J  JOuOOO'JOO 
-0.0o4o0b762 1 
-0. 1 1 74704766 
-0. 1585651^34 
-0. 18  /*52  7625 
-0. 2055733340 
-0.2114468578 
-0,205 5 7 33 34C 
-0.  lo7'>527625 
-0. 1585851434 
-0. L 174704766 
-0.06^608 7c2l 
-0. JOOOOOOOOO 
-0.0555555556 
-0.  Ill  I  1  11  11  1 
-0.151 7806004 
-0. 1 11  1  1  11  1  1  1 
-0.  1  HI  1  11  1  1  1 
-0.01<t68b0o62 


V2(  1) 

-o.ocooucoooo 

-o.ooocoooooo 

-o.oooooooooo 
-o.ocococoooo 

-0.0000000000 

-o.ccoccccooo 

-0.0000000000 

-o.cocoocoooo 

-0.0000000000 
-0.0380871927 
-0.0380871*2  7 
-0.  1523487/09 
-0.  I523<*d  //OS 
-C.  152348/709 
-  ).  152348770* 
-0.  15234e/709 
-0.  1523487709 
-0.152  3487709 
-0.03808  7192  7 
-O.C38C871927 
-0.0000000000 
-0.0  148860662 
-0.  lllll  1111  1 
-O.llllllUU 
-0.151  /806004 
-0.1111111111 
-0.0555555556 


buLUT  I  UN    U( X,Y ) 

0.  534722^222222 
0.256944 t4444444 
0.  2ul  3dacfc88888-> 
o. 36805555555556 
0.  4j055555555556 
0.  152 77777/77  778 
0. 097^^2222222 
0. 2&3LJb  Efc8888b* 
0.  ^40277  77/77778 
0. 0o250CC0000000 
0. 006*4444444444 
0.17361111111111 
0. 26;>6t,bcbd6ea8* 
-0.013888dBo8od89 
-  J. Jo944444444444 


0.44444444444444 
O.ZZ2ZZ2Z2ZZZZ<-t 
O.ZtZZZZZZZZZZZZ 
0.44444444444444 
0.340277/7///778 
0.  11805555555556 
0.1  1805555555556 
0.34027777777778 
0.25O0C0OCCC00O0 

0.02777777777778 

0.027/777  77/7778 

0. 250000000 COOOO 

0.  1  736111  11  lllll 

-0.0486111  1111111 

-0.04861111111111 


0. 36805555555556 
0.20138898888889 
0.25694444444444 
0.534  IZZZZZZZZZZ 
0.2636888868dd89 
0.097ZZZZ2ZZZZ2Z 
0. 15277777777778 
0.43055555555556 
0.  1736L111U1  11  1 
0.00694444444444 
0.06250000000000 
0.34027777777778 
Q.Q^IZZZZZZZZZZZ 
-0.06  9444  44444444 
-0.01388888888889 


0.30555555555556 
0.  19444444444444 
0.30555555555556 
0.63888888888889 
0.20138888888889 
0.09O27777777778 
0.20138888888889 
0.53472222222222 
0.11111111111111 

-0.00000000000000 
0.11111111111111 
0.44444444444444 
O.Oi^TZZZZZZZZZZ 

-0.07638888888889 
O.oi^lZZZZZZZZZZ 
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0.09122222222222  0.17361111111 

0.2013dbfc8886dd9  0  . 1 1 1  1 1 1 1  1 1  11 

-0.0/638do88888b9  -0. 1 1 1  1111 11 11 

-0.13194444444444  -J.  I  11  llll  11  11 

0.03^12222222^.^2  0.  1  1  1  1 1 1  1  1 1  1 1 

0.  152/77/////// 8  0.062 b 000  CO  CO 000 

-0.  12  5O0OCOO00000  -0  .1^12222222222 

-0.  Ib0b5b555555bb  -0  .lb9/  2222222222 

-0.01386888888889  0.06  25000000000  0 

O.lld05bbbbb55b6  0.027  7777/777/78 

-0.  159/2222222222  -0.  19444444444444 

-0.215277  77777778  -0. 194444444  44444 

-0.04  86  11  1  Llll  111  0.0277777//77778 

0.  09/ 222 2222 2 d^  0.006  9444 44 4444  4 

-0.  Id0bbbbb55  5556  -0.21b2  7/7  7  77777a 

-0.23  6  1  11  11  111  111  -0.21527777777  778 

-0. 06  94 4t 44444444  0.0069444  44  44444 

0.0  9  02777777/778  -0. 0000000 00 COOOO 

-0. 18 75 00 UC 000000  -0.222222222222*:2 

-0.2430 55 55555 556  -u  .22222222222222 

-U.0/b388fcfcd8  8bb9  -0 . 000000000 COOOO 

0.  09  I22c22 222 22^  0.006  9444  4  4  4444  4 

-0. Ia0  555555bbbbb  -0.215277/ 777777b 

-0.2361  I  1  I  1  111  I  11  -0.2152  777777  7778 

-0. 06  944444444444  0.0  06 ?444 4444444 

0.  UttOb5-3bbbbb:i6  J.  02777/777  Z/77o 

-0. 15972222222222  -0. 194  44444444444 

•0.21527777777778  -0  .  I  94  444444 44  4  44 

-0.04 doll  llll  1  I  1  1  0.02  777  7777 77ZZ8 

0.1b*://// 7  777/78  0  .O62bO000OCO0O  J 

-0. !2b0  0'JCu000000  -J.  159722222  22  222 

-O.ldO^bbbb^bbbSb  -0.  159/2  22  22  22  22*; 

-0.0lj>ada;:da3dBd9  O.O62500  0  00  00  000 

O.*:0l3  6dd£d6dB«9  O.llllllllllllll 

-0. 0/63  da  36838 H 39  -O.llllllllllllll 

-0.1jl94444  44-»4-4  -O.llllllllllllll 

0.U347222222  2222  O.llllllllllllll 

0.2o3bddedddaot>9  0.17361111111111 

-0.0l3HddSda838b9  -0.04d61lilllllll 

-0.06944444444444  -0.046b  1 1 1  1 1  1 1  I  I  I 

0.0W222<:2222222  0.173ellllllllll 

0.  J402  7  7  //  77  7  7  76  0.25000000000000 

0.062500      100  j.  J2 /////// Z/77o 

o. oob 944444 D. 02 777  777777778 

O.IM-  Llllllllll  0.25OJOOO0CCOOO0 

0.43  155555555556  ).34J277  77777778 

0.1^2  1  1 1  I  1  /t  1  lib  J. 11805555! 

d.0972Z2222Z2ZZ2  0.11805555555! 

0. 2b3odddda  0.   •    ' 1 1  I  t  1 1  1  lb 


0.2b3d8888888889 

0.03472222222222 

-0. 13  194444444444 

-0.0/638888888889 

0.20l38a88888889 

-0.01388888888889 

-0.18055555555556 

-0.  12500O000C0000 

0.152///77777/78 

-0.0486111  llll  HI 

-0.21527///////78 

-0. 15972222222222 

0.11805555555556 

-0.06944444444444 

-0.23611111111111 

-0. 18055555555556 

0.01122222222222 

-0.07638388888889 

-0.24305555555t>bo 

-0.18750000000000 

0.0902777//77778 

-0.06  94444  4444444 

-0.2361111111111 1 

-0. 16055555555556 

0.09/22222222222 

-0.048611 lllll 111 

-0.2152/77/  11111 6 

-0.  lb<5/2222222222 

0. 11805555555556 

-0.013ttBd8dd888d9 

-0.18055555555556 

-0. 1250O0OOOC0OO0 

0.15277777777778 

0.03472222222222 

-0.13  194444444444 

-0.0/638dde8888d9 

0.20  1388888888^  ^ 

O.09722222222222 

-0.0t>944444444444 

-0.01 338d38tfa8889 

0.2b3»868bd83889 

0. I  736111  11  1  I  II  I 

0.0069444444444-. 

O.Oo2bOOOOOOJOJ  J 

.  14 02  7777777778 

0.2  b 3383 68. 388389 

0.09122222222222 

>. 15277777777778 

.  .30555555555S6 


0.36805555555556 

-0.027777/////778 

-0. 13888888888889 

-0.027//7777////8 

0.30555555555556 

-0.0/638888888889 

-0. 18750000000000 

-0.07638888888889 

0.25694444444444 

-O.llllllllllllll 

-0.22222222222222 

-O.llllllllllllll 

0.22222222222222 

-0. 13194444444444 

-0.24305555555556 

-0.  13194444444444 

0.20138888888889 

-0.13888888888889 

-O.250OOO00C00O00 

-0.  13888886833889 

0.  19444444444444 

-0.  13194444444444 

-0.24305555555556 

-0.  13194444444444 

0.20138888888889 

-O.llllllllllllll 

-0.22222222222222 

-O.llllllllllllll 

0.22222222222222 

-0.0/638888833889 

-0.18/50000000000 

-0.0/638838888889 

0.25694444444444 

-0.O2///////7///8 

-0.  I  3688888d8dd89 

-0.02///////////8 

0. 30555555555556 

0.03^12222222222 

-0.0/638688888889 

0.034Z2222222222 

0.  36805555555556 

O.llllllllllllll 

-0.00000000000000 

O.llllllllllllll 

0.4444  •  •  ^4444444 

0.201 3883888a639 

0.0902// 7  777/7/8 

0.201 3d888886dd9 

0.  '33^12222222222 


13  8       -I.OOOOluuolO  1.73^0506076         0.0 

13  9       -0.9160797631       -2.3166247904  0.2000000000 

13  10       -0.656b542495       -1. 316624 79C4  0.4CO0C0OO00 

13  11       -0.1961524227       -0. 3 16624 7904  0.6000000000 

12  12       -0.472135955C       -0. 4 72 1 3595 5C  l.OOOOCOOOOO 

11  13       -0.3166247904      -0.196152  422  7  C.6000C00000 

10    13   -l.316fc^47s04   -0.6568542495    0.4GCC  ^0000 

9  13       -2.3166247904       -0.  9160 79783 1  0.2  000 COOoC "! 

8  13  1.3024355C95      - 1 . COOOCOCOCC  0.0 

7  13  0.3024355095       -U.  432o  1 1  ol  7  7         0.7710628438 

7  12  0.7562036626       -  I . 4328  1 1 fl 1 7 7  C. 6168502751 

6  11  0.040524284b       -0.1720  304203  0.2313188532 

6  10  0.2214414fc91       -1.1720  304203  0.1542125688 

6     9    0.3232390981   -2.1720304203    0.0771062844 

6  8  0.3561944902  3.1720304203  0.0 

6  7  0.3232390981  2.1720304203         0.0771062844 

6  6  0.2214414691  1.1720304203         0.1542125688 

6  5  0.0405242648  0.1720304203  0.2313168532 

7  4  0.7562036828  1.4328116177  0.6168502751 

7  5  0.3024355095  0.4328118177  0.7710628438 

8  t>       -1.  7320506C76  l.OGOOOOCOOO         0.0 

9  3       -0.7320508C76         0.7220508076  1.0000000000 
9            4       -l.COOOOOOCOO          1.7320506076         0.0 

10  5       -0.2679491924  1.0000000000  0.5000COOOOO 

11  6       -l.OOOOOOOOOU  0.^679^91924  C.50OOC00000 

12  7       -1. 7320508G76  1.0000000000         0.0 

13  7       -0. 732C506C76         0.732C506076  l.OOOOCOOOOO 
0.0                                  0.0 

0.0  0.0 

0.0  0.0 

0.0  0.0 

0.0  0.0 

0.0  o.c 

0.0  0.0 

0.0  0.0 

0.0  0.0 

-0.0646067621  -0.03808719^7 
■0.1174704766  -0.03  808  7192  7 
-0.1585851434  -0.152348  7709 
-0.1879527625  -  0.  1  5234t>7709 
-0.2j55733j4C  -  0.  1  52 348 77C9 
-0.2114468578  -0.152348  7709 
-0.2055733340  -0.1523467709 
-0.18  79527625  - C. 1 52348 7709 
-0. 1D85851434  -0.1523487709 
-0.1174704766  -0.0380871927 
-0.0b460676^!l       -0.0380871927 

0.0  0.0 

-0.0555555556       -0.0148  860662 

-o.iiiiiiiui     -o.il  u  linn 

-0.1517806004  -0.  1  1  1 1 1 1 1 1  1 1 

-0.  11  11  1  lllll  -0.  1517606004 

-O.lillllllll  -0.1111111111 

-0.0143b60o62  -0.0555555556 

0  0.0  2.62C52D+CI     1.459700*00 

1  -4.4657021980D-01  1.814980*00  <*.65756D-01 

2  -6.2448872675D-C1  9. 117440-01  2. 216750-01 

3  -5.64799223230-01  5.00C49D-01  1.203670-01 

4  -6.4129364256U-C1  3. 750620-01  6.700450-02 

5  -6.0671  lb9046D-01  4.093940-01  2.826410-02 
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6  -5.89888503680-01  2.  1027GD-01  1.C8533D-02 

7  -5.6516813624D-C1  9.043260-02  2.904030-^ 

8  -5.41701078710-01  2.951120-02  6.176290-04 

9  -5.35067680820-01  5.96C87D-03  1.70681D-04 

10  -5.33920436550-01  1.962850-03  3.258720-05 

11  -5.3437389326D-C1  4.816380-04  7. 683310-06 

12  -5.3453890774D-01  1.254650-04  1.546630-06 

13  -5.3457758425C-01  2.21249D-05  3.196700-07 


NO.  OF  ITEkATlLNS: 

STOPPING  CONDITION  *AS: 

NOKM  OF  RESIDUAL  «AS  REQUIRED  TO  BE 

ESTIMATE  OF  REASONABLE  NORM: 

ESTIMATE  OF  NORM  ACTUALLY  OBTAINED: 


13 

1 
1.00000D-06 
2.35465D-13 
3.  1V670D-O7 


ESTIMATES    OF     NORM    OF     FINAL    RESIDUAL 

LQ:  3.19670D-07 

LO  INCOMPLETE:     2.21249D-05 

QK:     1.514510-06 
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